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ABSTRACT 

The  basic  outline  of  a  digital  computer  (CRC-102A)  program  for  the 
numerical  solution  of  the  approximation  problem  in  the  time  domain  is 
presented.   The  method  used,  in  effect,  is  a  real  time  convolution  and 
Dirichlet  series  representation  of  the  Laplace  transform  to  effect  a 
time-to-frequency  transformation.   The  major  portions  of  the  digital 
computer  program  and  the  complete  flow  chart  to  solve  the  problem  are 
included. 

The  composition  of  a  basic  library  of  programs  for  the  solution  of 
all  network  synthesis  problems  and  the  integration  of  this  program  into 
such  a  library  is  discussed. 
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1.   Introduction 

The  problem  of  network  synthesis  is  a  difficult  and  complex  one. 
Consider  the  synthesis  process  to  be  divided  into  the  two  separate  sub- 
problems  of  approximation  and  realization.   There  are  many  techniques 
for  the  attainment  of  the  first  goal,  Z(p),  the  transfer  function. 
Similarly,  having  attained  Z(p),  there  are  many  procedures  available  to 
realize  Z(p)  in  any  one  of  several  different  circuit  configurations. 
From  the  point  of  view  of  the  electronics  research  engineer,  a  tool  which 
would  rapidly  and  accurately  compute  the  circuit  components  of  a  network 
to  provide  a  specified  response,  at  a  specified  impedance  level,  in  a 
specified  circuit  configuration  would  free  the  researcher  from  the  tedium 
of  designing  the  network  and  also  from  the  temptation  of  accepting  a 
"second-best"  solution.   The  designer  would  feel  no  compulsion  to  term- 
inate a  problem  rather  than  face  another  "messy"  synthesis  problem. 
The  modern  high-speed  digital  computer  offers  the  engineer  such  a  design 
tool,  if  an  appropriate  library  of  synthesis  programs  is  available. 

Possibly  the  greatest  advantage  the  computer  avails  its  users  is 
that  once  a  certain  type  of  problem  has  been  solved,  the  method  of  solu- 
tion may  be  permanently  recorded  on  punched  cards  or  on  magnetic  tape  or 
paper  tape.   The  engineer  now  need  only  recall  that  such  a  problem  has 
been  programmed  and  refer  to  the  catalogue  of  computer  programs  rather 
than  to  a  textbook  to  refresh  his  mind  on  how  to  solve  the  problem  by 
paper  and  pencil. 

For  a  comprehensive  summary  of  approximation  methods,  the  reader 
is  referred  to  Trans,  IRE,  Prof  Grp  on  Circuit  Theory,  Vol.  CT-1,  No.  3, 
of  Sep  1954.   In  addition  to  the  summary,  the  first  paper  presented  by 
Stanley  Winkler,  has  a  definitive  bibliography  of  240  items  on  network 
synthesis. 


It  is  the  purpose  of  this  paper  to  present  a  program  for  the  solu- 
tion of  the  approximation  problem  in  the  time  domain  and  to  show  how  this 
program  may  be  integrated  into  a  library  of  programs  which  use  convention- 
al methods  for  the  solution  of  the  more  common  network  synthesis  problems. 


2.   Statement  of  the  Problem 

The  approximation  problem  as  considered  in  this  paper  is  felt  to  be 
defined  as: 

Given  the  excitation  available  or  expected  at  the  input  terminals 
of  a  passive  four-terminal  network,  and  the  response  or  output  required, 
or  the  Impulse  response,  where  these  parameters  are  expressed  in  graphical, 
or  general  analytic  form;  find  a  rational  function  expressed  as  two  poly- 
nominals,  where  one  is  understood  to  be  the  numerator  and  the  other  the 
denominator,  which  describes  the  transfer  function  within  the  limitations 
on  the  allowable  deviation  of  the  output  or  on  the  number  of  elements 
to  be  used  (but  not  both),  specified.   Provisions  are  to  be  included  for 
compensation  of  dissipatlve  effects.   The  polynomlnals  will  be  expressed 
either  as  polynomlnals  or  as  the  quadratic  factors  of  polynomlnals.   The 
function  is  to  be  physically  realizable. 


3.   The  Approximation  Problem  in  the  Time  Domain 

Heretofore,  approximation  problems  in  the  time  domain  have  involved 
lengthy  and  inaccurate  frequency  domain  expansions  which  do  not  always 

lead  to  realizable  system  functions.   These  methods  involve  complicated 

1(3) 
integrations,  or  at  best,  Taylor  series  expansions.      At  any  rate, 

they  are  not  suitable  for  the  flexible  program  requirements  we  desire. 

What  we  would  like  to  have  is  a  method  which  would  rapidly  transform  the 

input-output  characteristics  of  one  network  into  a  realizable  system 

function  within  any  desired  degree  of  accuracy,  or  alternatively,  an 

optimum  realization  for  a  specified  number  of  elements.  The  numerical 

1(3) 
method  of  Ba  Hli     is  tailor-made  for  these  requirements.   Furthermore, 

not  only  will  it  operate  on  input -output  relationships,  but  also  there 
is  no  reason  why  an  impulse  response  which  has  been  obtained  analyt- 
ically cannot  be  plotted  and  introduced  into  the  program  at  the  appropri- 
ate point  to  produce  the  transfer  function  for  all  types  of  networks, 
i.e.,  predictor  networks,  smoothing  filters,  interpolation  filters,  com- 
pensation networks  and  so  on.   In  fact,  we  may  say  if  a  physically  realiz- 
able h(t)  can  be  defined,  we  can  find  the  system  function.  The  designer 
must  have  a  plot  or  sequence  of  ordinates  which  describe  the  h(t),  or 
the  f,(t)  and  f  (t),  and  use  this  plot  or  sequence  to  enter  his  problem 
into  the  approximation  program. 

Briefly  stated,  Ba  Hli's  method  is  to  find,  by  a  process  of  synthetic 
division,  a  sequence  (  \hj     )  which  represents  the  impulse  response  of 
the  network  which  in  turn  relates  the  input  waveform  to  the  desired  output; 

Freddy  Ba  Hli,  Network  Synthesis  by  Impulse  Response  for  Specified 
Input  and  Output  in  the  Time  Domain,  Tech.  Rpt.  No.  261,  MIT  Res.  Lab. 
of  Electronics,  July  31,  1953  (hereinafter  abbreviated  as  NSIR-BH). 


and  then  to  calculate  the  psuedo -system  function  H*(s),  from  Vhj   by  a 

simple  power  series  expansion.   The  H*(s)  is  a  polynominal  in  s  which  is 

Pm(s) 
the  ratio  of  the  Qn*  f  which  is  normally  considered  to  be  the  expression 

of  the  system  function.   In  effect,  we  convert  the  convolution  integral 

into  its  component  operations  and  use  an  iterative  substitution  method. 

We  shall  not  reproduce  the  theoretical  reasoning  and  proof  which  was  so 

neatly  done  by  Ba  Hli,  but  we  shall  demonstrate  the  ease  and  simplicity 

of  the  method  with  a  numerical  example. 

Before  proceeding  with  a  complete  example ,  let  us  show  how  the 

synthetic  division  process  yields  the  impulse  response  in  a  couple  of 

trivial  cases,  the  results  of  which  are  well-known.   Consider  Figure  1, 


.5 
0 


Unit  Step 
(f^t)]   -  1,  1,  1,  1,  1, 


Unit  Ramp 
[fo(t)j   -  .2,  .4,  .6,  .8,  1,0, 


A  Trivial  Case 
Figure  1 


the  input  f.(t)  we  take  as  the  unit  step  whose  Laplace  transform  is  — j 
1  s 

£  (t),  the  output,  is  the  unit  ramp.   Its  transform,  of  course,  is  —2. 
o  s 

If  we  sample  the  input  and  output  at  intervals  of  .2  of  a  time  unit,  the 

sampled  ordinates  for  f.(t)  are  A  1, 1,1,1, 1,1, ... ,j,  and  for  f  (t) 

(.2,  .4,  .6,  .8, 1.0, 1.4, 1.6 "I.   Now  dividing  (f  (tj)/      in  the  manner 

,  ,  i^i  ta  i  i^i |2,  «2lt.ttt 


illustrated:   1, 1,1,1,1, 1. .../. 2, .4, .6, .8, 1.0,1.2,1.4 

we  have  \hjA  ■  .2,  .2,  .2,  .2,  .2 and  if  we  divide  \h\  by  At,  we  pro- 
duce the  sequence  1,1, 1, 1,1, ... .again  the  unit  step.  This,  of  course,  is 
completely  analogous  to  the  process  in  Laplace  transform  operations,  viz: 

or     H(s).Fo<^(8)       .10i.   .1/. 

We  can  apply  the  process  in  reverse,  for  example,  jf  (t)r  x  ^hj 
■  jf  (t)l,  if  we  use  the  unit  doublet  as  \hj    we  should  differentiate 
the  output.   Consider  the  sequence  below  which  defines  a  sine  wave,  with 
At  -  .1 

{f^  -  .1   .199   .296   .389   .479   .564   .644   .717 
[hV  -  10    -10 


1    1.99   2.96   3.89   4.79   5.64   6.44   7.17 

-1     -1.99   -2.96  -3.89   -4.79   5.64   6.44 

[fjj    1     .99    .97    .93    .90    .83    .80    .73 

[fj  (Cont'd)   .783   .841    .891   .932   .964   .985   .997 

7.83   8.41   8.91   9.32   9.64   9.85   9.97 
7.17   7.83   8.41   8.91   9.32   9.64   9.85 

}fQ]  (Cont'd)    .66    .58    .50    .41    .32    .21    .12 

6 


and  if  (t)f  is  a  sequence  which  defines  the  cosine  which  is  as  expected 
since  the  unit  doublet  acts  as  a  differentiator.   Thus  we  demonstrate 
the  ability  of  this  S^bi    to  perform  as  an  operator  equivalent  to  H(s) 
in  the  frequency  domain  on  functions  defined  as  sequences  in  the  time 
domain.   The  transformation  of  {h)   to  H*(s)  follows  logically  from 

the  definition  of  the  convolution  integral: 

t 
fQ(t)  -  J   £t(t)  h  (t  -T)  dT 

(3) 

Ba  Hli     has  shown  that  H*(s)  can  be  represented  by  the  Dirichlet 

power  series, 

H*(s)  -  J       a  exp  (-Tns) 
n-1 

where  the  a  are  the  areas  in  our  |_hj   sequence,  and  the  T  are  de- 


fined by 

h^     k  .    A  t 

I   "  n/\t r— 

n  2   , 

that  is,  the T  are  measured  from  the  origin  to  the  center  of  the  sampl- 
ing interval  under  consideration  at  the  moment. 

We  may  approximate  H*(s)  as  closely  as  we  please  by  taking  our 
sampling  points  at  smaller  and  smaller  intervals. 

To  demonstrate  the  complete  process  and  provide  a  basis  for  discuss- 
ing some  of  the  more  obscure  points,  let  us  take  another  example  and 
follow  it  through  more  closely.   Figure  2  shows  the  f  (t),  the  desired 
f  (t)  and  the  resulting  h(t)  required,  (At  «  0.1): 
{^(t)}-  1,    1,    1,    1,    1,   0,   0,   0, 
{fo(tjj-.l,   .2,   .3,   .4,   .5,   .4,   .3,   .1,   0,   0,   0, 
\h]  A  -.1,   .1,   .1,   .1,   .1,   0,   0,   0,   0, 
In  tabular  form,  we  have 


Tn 
.05 
.15 
.25 
.35 
.45 
Substituting  Into  our  expression  for  H*(s),  we  have 


H*(s) 


[l  -  .05 

\}  '  -15 

[l  -  .25 

[l  -  .35 

[l  -  .45 


s  + 


s  + 


s    +• 


s    + 


s    + 


(.05s)4 
2! 

(•15s): 
2! 


2! 


(r35s)' 
2! 


3! 

C15s): 
3! 


(«25s)- 
3? 


31 


(.45s)- 
2! 


U£*X 

3! 


(.05s) 

4! 

4! 


(r25s) 

4! 


(r35s) 

4! 

U££L 

4.' 


-■■•] 
...3 

-••3 


-   .5  -  .125s  +  .020625s   -  .00255208s  +  ,0002517968s   -  

If  we  compare  this  result  with  that  given  by  Laplace  transform  methods, 
we  have 

^[fi(t)]   ■  ~  [l  -  exp  (-l/2s)] 

^[f0(t)l   -  72  C1  "  exP  (-l/2s)]  2 

2[h(t)]   ^&-(t^|?l(^-i[l-«p(.i/l.l 

and  the  power  series  expansion  for  ^\h(t)J   is 

.5  -  .125s  +  .020833s2  -  .002604166s3  +  .0002605166s4 
which  shows  the  first  two  terms  by  Ba  Hli's  method  agree  exactly,  a  IX 
error  in  the  third  term,  2%  in  the  fourth,  and  about  47.  in  the  fifth. 


Before  proceeding  let  us  see  how  this  accuracy  may  be  improved. 
Clearly  increasing  the  number  of  terras  taken  and  maintaining  the  same 
interval  will  have  no  effect  on  the  formation  of  the  second,  third,  and 
fourth  degree  terras  whose  accuracy  we  wish  to  improve  and  intutively  one 
feels  that  increasing  the  number  of  sample  points  will;  therefore  let  us 
increase  the  number  of  sample  points  to,  say,  8. 

With  f.(t),  f  (t)  and  h(t)  as  before  we  tabulate,  for  At  -  .0625: 


Tn 

.03125 
.09375 
.15625 
.21875 
.28125 
.34375 
.40625 
.46875 


n 

.0625 
.0625 
.0625 
.0625 
.0625 
.0625 
.0625 
.0625 


H(s)  -   .0625     [l   -   .03125s  ♦   <'°>™*>     -  UfllpJL  +  (-03125,)     .    _J 

.0625  [l  -  .09375s  +  (-09375.)2  .  (,09375s)3  .  (-09375s)4  .  _J 
and  so  on;  the  resulting  expression  is: 

H*(s)  -  .5  -  ,125s  ♦  .020635719s2  -  .002631105s3  +  .000257032861s4  -  ... 
and  the  new  expansion  for  8  sample  points  has  produced  errors  of  0%  in 
the  first  two  terras,  .95%  in  the  third  term,  1%  in  the  fourth  and  1.3% 
in  the  fifth.   Therefore,  the  increase  in  the  number  of  points  has 
obviously  resulted  in  an  overall  improvement.   A  further  increase  to 
10  points,  produces  accuracies  of  0%,  07.,  .25%,  .5%,  .8%  and  1.25%. 

The  foregoing  results  are  tabulated  in  Table  1  for  the  purpose  of 

10 


easy  comparison. 

We  recognize  that  the  errors  in  the  higher  order  terms  are  caused 
by  replacing  the  Laplace-Stieljes  integral   by  a  truncated  series  summa- 
tion.  This  is  equivalent  to  saying  that  we  have  not  computed  the  true 

'•moments"  of  a  about  the  origin, 
n 

Consider  Figure  3,  we  have  here  a  triangular  impulse  response.   Let 
us  calculate  the  second  and  third  moments  of  this  curve  and  compare  the 
results  with  the  power  series  expansion  of  the  Laplace  transform  and 
also  the  Dirichlet  series  expansion.   We  let  b_  be  the  second  "moment", 

(it  would  be  more  precise  to  call  this  by  a  different  term  since  it  is 

M 
actually   1  .  where  M   is  the  moment),  b   equals  the  third;  thus, 

j!        j  3 


1.0 


h(t)     .5 


h(t)  -  t,      0  ^  t  ^  1/2 
1-t,     l/2^t^l 
0,      elsewhere 


A  Triangular  Impulse  Response 
Figure  3 


^SIR-BH^,  pg.  15. 
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|y    L    0J  xdx+      1/2f         (1-x)   x  dxj 


1/2 


"     2!     L-4    >0  3      1/2       A    'l/2J 


1  i_     l     i_     I     L_v 

2  C64  +  3   "  24   "  4  +  64; 


.03645833 

*   r    /1/2    *a 

IT   L     0j  x  dx 

ci  1/2 


+     l/2f        (1_x)  x  dxJ 


3!     L5    '0  4     1/2  5    'l/2J 

1  JL_  I  1_  1  1       V 

"     6      U60  +  4   "  64   "  5  +   160' 
-     .00781250 

The  trancendental  Laplace  transform  for  this  function  Is: 

H(s)  -j[l-  exp  (-l/2s22 
s 

whose  power  series  expansion  is 

H(s)  -  .250  -  0.1250s  +  .03645833s2  -  .00781250s3 

The  Dirichlet  series  results  from 


•n 

n 

.1 

.02 

.3 

.06 

.5 

.09 

.7 

.06 

.9 

.02 

*TnS 
H(s)   -  y    an  e 

n-1 

Therefore,  considering  only  the  b  and  b.  terms,  we  compute 
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2  2  2  2  2 

b2  -  .02  x  ■—-  +  .06  x  -yr  +  .09  x  -yr  +  .06  x  yy  ♦  .02  x  yy 


j  (.0002  +  .0054  +  .0225  +  .0294  +  .0162) 


.03685 


Similarly,  b.,  is  computed  as 


b  »  .008008 
Therefore: 


••Moment" 

.03645833 

.00781250 


Ba  Hli  Approx. 
.036850... 
.008008... 


Laplace 
.03645830 

.00781250 

j 

The  obvious  method  of  reducing  the  error  is  to  take  a  smaller  in- 
terval and  thereby  approach  the  true  integration  to  a  closer  degree  and 
since  this  will  be  a  one  time  only  calculation,  we  are  not  adamant  about 
making  extra  computations  to  improve  the  accuracy.   We  also  notice  that 
by  making  our  intervals  smaller  we  reduce  those  errors  which  may  arise 
from  the  process  in  which  we  obtain  the  \h\      .   For  a  more  complete  dis- 
cussion of  the  reduction  of  the  inherent  error  in  Ba  Hli's  method  see 
Appendix  II. 

Having  obtained  a  psuedo-system  function,  it  yet  remains  to  deter- 
mine its  equivalent  as  the  quotient  of  two  polynominals,  the  degree  of 
each  as  yet  unknown.   We  equate  our  psuedo-function,  H*(s)  to  an  express- 
ion of  the  form 


H*(s) 


P£s2 
Q(s) 


pQ  +  Pls  +   ...    pms 


m 


n-1 


+  q   s 
n 


qQ   +  qLs   +    ...    q^s 

and  solve  for  the  unknown  p's  and  q's,  since  these  coefficients  express 
a  ratio  we  may  arbitrarily  set  any  one  of  them  we  choose,  (normally  q  ) 
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equal  to  one  without  any  loss  of  generality. 

We  notice  that  the  series  of  equations  which  results  from  this  opera' 
tion  will  always  have  the  form: 

po  a  hoqo 

pi  -  hiqo  +  Vl 

p2  -  h2qQ  +  h^  +  hQq2 

0   -  h3q0  +  h2qL  +  hLq2  +  h^ 

0   -  b4qQ  +  h3qL  +  h2q2  +  h^   +   h^ 

0   -  h5qQ  +  h4qL  +  h3q2  +  h2q3  +  h^ 

0   «  h6qQ  +  h5qi  +  h4q2  +  h3q3   +   h2q4 

For  m  ■  2,  n  *  4;  r,  the  relative  degree  between  P(s)  and  Q(s),  is 
obviously  two.   We  see  that  we  have  seven  equations  in  eight  unknowns, 
which  is  as  it  should  be  for  the  expression  of  a  ratio,  and  hence,  re* 
quiring  that  we  set  q.  (or  any  other  coefficient)  equal  to  some  conveni- 
ent value.   The  solution  of  this  portion  of  the  problem  is  considered  in 
detail  in  Appendix  III. 

We  choose  r,  the  relative  degree  of  numerator  and  denominator  by 
observing  the  behavior  of  h(t)  around  0.   Since  we  know  from  the  initial 
value  theorem  of  Laplace  Transform  Theory, 

lim   h(t)   -  sH(s) 

t — ?0         s — =>^ 

we  may  say  that  the  function  must  behave  as 

m 
p  s 

as  t  approaches  zero. 

vm 

M.  F.  Gardner  and  J.  L.  Barnes    ,  page  265. 
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Thus,  If  h(t)  starts  out  as  a  step 

P(s)  / 
v  /ft,    v — *   1/s  as  s  approaches  infinity;  or  is  h(t)  starts 

out  as  a  ramp       /  c\t    \ — *  ^s  as  s  aPProacnes  infinity. 

To  program  the  computer  to  sense  this  limit  and  thus  decide  whether  h(t) 
is  an  impulse,  step,  ramp,  or  higher  order  function,  we  first  determine 
h  ,  i.e.,  h(0).   Since  in  some  cases  this  may  result  in  division  by  zero, 
we  perform  a  first-order  backward  interpolation,  using  h  and  h  to  find 
h..   This  is  approximate  but  certainly  is  sufficiently  accurate  to  deter- 
mine whether  or  not  h(t)  starts  as  a  Jump  function  or  as  a  ramp,  which  is 
the  primary  decision  in  the  "selection  of  relative  degree  routine"  (see 
Appendix  I,  page  1-1).   If  it  is  decided  at  this  point  that  h(t)  is  a 

first -order  approximation  to  a  ramp  function,  we  must  further  decide 

3     4 
whether  it  is  truly  a  ramp  or  a  higher  order  function,  i.e.,  1/s  ,  1/s  , 

etc.   These  decisions  are  based  on  the  fact  that  the  n —  derivative  of  an 
n-order  polynominal  is  a  constant,  the  process  is  flow  charted  for 
+1  ^.  r  ^  -4  on  page  1-2-1  of  Appendix  I.   The  routine  may  be  easily  ex- 
tended to  encompass  all  possibilities  which  may  arise. 

We  know  that  we  can  obtain  any  desired  degree  of  accuracy  in  the 
terms  of  our  power  series  expansion  and  further  that  we  can  correctly 
determine  the  relative  degree  between  our  numerator  and  denominator  poly- 
nomials.  We  have  shown  that  we  can  compute  the  coefficients  of  our  H*(s) 
and  we  know  that  since  this  is  a  ratio  with  arbitrary  coefficients  we  may 
multiply  the  expression  by  the  impedance  level  factor  to  set  any  required 
impedance  level.   But  what  of  the  approximation  error,  i.e.,  the  error 
caused  by  not  taking  an  infinite  number  of  terms  to  approximate  tran- 
cendental  functions? 

Ba  Hli  has  shown  that  in  the  case  where  the  Laplace  transform  of 
h(t)  is  rational,  it  will  be  recovered  exactly.   NSIR-BH,(3)  pgs.  41-44. 
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N  »  5  -1  +  1.5(5)  +2   -   13.5,  which  we  take  as  14. 
Therefore,  we  can  write 

P0  +  P1S  +  P2g2  *   •••  *   P12s12  ^P13s13    u         .  .    11       K   27 

2 13 14  "  h0  +  hlS  +  *••  hllS   +  ••'  h27S 

q0  +  ql8  +  q2S  +  •••  +  qi3s   +  qi4s 

this  expression  leads  to  23  equations  in  24  unknowns,  one  of  which,  custo- 
marily q~,  we  set  equal  to  1,  without  error.   To  avoid  filling  sheets  of 
paper  with  useless  symbology,  we  shall  terminate  this  example  here  and 
turn  to  another  pair  of  functions. 

We  shall  compute,  using  the  program,  approximate  impedance  functions 
for  a  pair  of  input -output  time  functions  and  check  the  computer  solution 
analytically.   To  do  this,  we  choose  the  functions  illustrated  in  Figure 
4. 

f.(t)  is  a  rectangular  pulse 
f  <t)  -  1    0  £  t  ±  1 
0    elsewhere 

f  (t)  is  an  inverted  sectioned  parabola 

f  (t)  -   .5t2  0  ±   t  **   1 

o 

-  -t2  +  2.0t  -  1.5     1  ±   t  ±  2 

-  .5    (t-3)2  2it^3 
a     0  elsewhere 

for    At   ■    .1,    the   following   sequences   obtain 

^(t)]      -     1,  1,  1,  1,  1,  1,  1,  1,  1,  1, 

{fo(t)j      -   .005,    .020,    .045,    .080,    .125,    .180,  .245,  .320,  .405,  .500, 

.590,    .660,    .710,    .740,    .750,    .740,  .710,  .660,  .590,  .500, 

.405,    .320,    .245,    .180,    .125,    .080,  .045,  .020,  .005,  .000, 

Synthetic  division  of  jf   (t)/f   (t)l  yields : 
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[h]    -  .005,  .015,  .025,  .035,  .045,  .055,  .065,  .075,  .085,  .095, 
.095,  .085,  .075,  .065,  .055,  .045,  .035,  .025,  .015,  .005 
which  describes  h(t)  as  the  triangular  function  also  shown  in  Figure  4. 
We  know  that  the  Laplace  transform  from  f.(t)  is: 

^[fi(t)]    "     8    (1  '  C"S) 

and  the  Laplace  transform  for  f  (t)  is: 


-s-,3 

-S-.2 


Thus  H(s),   which  should  be     [*   <t)]      /      [j^CO]    °r  [*   "  e    j    ,   which  is 

2 
the  transform  of  a  triangular  pulse  . 

We  expect  H*(s)  as  computed  by  our  program  to  be  a  good  approximation 

n  -  e"81  2 
to  the  power  series  expansion  of  

r^p — J   -  1.000000000   -  1.000000000s   +  .583333333s2 

-  .250000000s3  -   .086111111s4   -  .025000000s5 

+  .006299603s6  -   .001405423s7  +  .  .  .  . 

The  program  approximation  with  single  interpolation,  i.e.,  At/2,  produced: 

H*(s)  »  .999999999    -  1.000000975s   -  .583340832s2  -  .250007019s3  + 

.086115077s4  -   .025001670s5  -  .006300172s6  -  .001405587s7  +  . 

If  we  now  solve  for  the  impedance  function  using  n  »  3,  and  r  =»  -2 

(since  h(t)  is  a  ramp  function  in  the  vicinity  of  zero); 

,,  x     0.999999998  -  0.049971105s 

Z(p)  »  -* 2 T    or; 

1.0  +  0.950029868s  +  0.366689961s  +  0.062506119s 

-     (-.049971105)   (s   -  20.011564640) 

(s  +  2.277182620)      [(s  +  1.79464120)2  +  1.9505899912] 

Gardner  and  Barnes    ,  Appendix  A,  No.  6.08. 

2  (  k  \ 

Gardner  and  Barnes    ,  Appendix  A,  No.  6.07. 
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The  inverse  Laplace  transform  of  this  function  is: 

4.4132006227  e"2»27718262t  +  4.465570186  e-l-79464120t  sln  (l.95058999t-81.217°) 
The  plot  of  this  inverse  superimposed  on  the  desired  h(t)  is  shown  in  Figure 
5a. 

If  we  desire  to  improve  this  approximation  we  go  to  a  4-pole  network. 
The  program  computes  the  system  function  as: 

2.67414555s2   -  8.1622705s  +   135.9105227 

s4  +   10.59828491s3  +  51.1404969s2  +  127.7483847s  +  135.9105230 

which  we  expand  as : 

1.0681279249  {- <9  »  23,98432161)  .     (s  +  20.00984647) \ 

I  [(s+3. 01 7294581)  +(.932721456)  J         [( s+2. 28184 7873)  +(2.901655110)  JJ 

and  determine  the  inverse  Laplace  transform  to  be: 

24.010971222  e*3-017294581t     sln   (.932721456t  +  2.547°) 

-6.612687409  e-2»281847873t   8in  (2.901655110t+  9.295°) 
This  inverse  is  plotted  in  Figure  5b.   Note  the  improvement  with  the 

addition  of  the  other  pole.   The  approximation  can,  of  course,  be  further 

improved  by  adding  still  more  poles. 

Suffice  it  here  to  say,  that  the  existing  program  will  solve  for  the 

impedance  function,  given  the  time-series  data  describing  the  input  and 

output  waveforms. 

Before  concluding  this  section,  we  estimate  the  required  solution 

time  for  the  total  approximation  problem  in  the  time  domain. 

We  assume  the  following  conditions  obtain  as  a  basis  for  our  estimate: 
1.   h(t)  is  not  a  closed  function,  i.e.,  K  ■  4  (see  Appendix  II).   If 

h(t)  is  closed,  that  is,  if  h(t)  is  completely  described  in  4  or  less  time 

units  the  running  time  for  the  synthetic  division  routine  and  for  the  com- 
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putation  of  the  psuedo-system  function,  H*(s),  will  be  proportionately 
reduced. 

2.   H(s)  starts  as  a  step  function,  i.e.,  r  -  -1.   This  is  the  worst 
case  we  encounter  as  regards  computation  time  for  the  psuedo-system  funcf 
ion,  since  we  must  compute  2(n  +  1)  +  r  terms  of  the  expansion  of  H*(s) 
in  order  to  solve  the  system  of  simultaneous  linear  equations. 

3.   The  interpolation  decision  is  such  that  one  intermediate  inter- 
polation calculation  is  performed.   This  is  the  intermediate  condition, 
and  requires  twice  the  computation  time  for  a  no-interpolation  decision 
and  one-half  the  time  for  the  four-point  interpolation. 

We  consider  the  overall  program  (for  time  domain  approximation)  to 
be  sub-divided  into  three,  (four  if  Z(p)  is  desired  in  factored  form  or 
if  compensation  is  to  be  provided)  major  computations.   Thus  we  have. 
Sub-division  CRC  time  Word -times 

Synthetic  Division  to        5  min*  72  x  in4* 

produce  (hi  "■   x  i0 

Dirichlet  power  series       n  x  3  min  (n-D  x  44  x  lfl4 

expansion  H*(s)  l    '  x  **   x  i0 

Matrix  Inversion  to  .064  n3  min  96  n3  x  104 

compute  Z'p  *  xu 

Factor  Z'p  to  apply  (2n-l)  x  9.4  min        (2n-l)  x  144  x  104 

loss  compensation  v    }  H   x  i0 

♦Estimated,  all  other  values  measured. 

For  a  10-pole  network,  this  gives  a  solution  time  of  278  minutes;  for 
a  15-pole,  538  minutes;  for  a  20-pole,  945  minutes. 
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4.   The  Network  Synthesis  Library 

Figure  6  shows  the  block  diagram  of  a  network  synthesis  library  for 
the  automatic  computation  of  the  approximation  functions  for  a  variety  of 
filters.   This  diagram  is  by  no  means  considered  to  be  all-inclusive,  it 
is  only  intended  to  show  how  such  a  library  would  be  assembled. 

It  should  be  obvious  from  the  diagram  how  our  various  blocks  are  used 
together  to  provide  a  large  number  of  filter  characteristics  from  a  rela- 
tively small  number  of  programs.   All  blocks  would  have  common  input, 
output  and  carry*  locations,  and  programs  represented  by  the  blocks  would 
have  bootstrapping  links  to  call  in  the  next  block  in  order  to  proceed 
with  the  computation. 

In  addition  to  the  "frequency"  filters  shown,  we  would  consider  it 
desirable  to  provide  in  the  basic  library: 

a.  Tschebycheff  Stop-Band  Ripple  Filters 

b.  Tschebycheff  Equi-Ripple  Filters 

With  regard  to  the  latter,  Byron  J.  Bennett   has  developed  an  inter- 
esting iterative  technique  which  seems  to  be  much  simpler  than  Alexander 

2 
J.  Grosman's  method  for  avoiding  the  elliptic  function. 

Note  the  provision  for  the  realizibility  check.   It  is  considered  that 

the  inclusion  of  such  a  routine  is  mandatory  in  order  to  save  computer 

3 
time.   This  routine  would  compute  n,   the  number  of  elements  required  to 


*We  use  the  term  carry  to  characterize  that  data  which  is  not  used 
in  the  present  computation  but  must  be  "carried"  for  use  in  later  routines 

Byron  J.  Bennett,  A  Note  on  Circuit  Synthesis,  p.  61. 

2 
A.  J.  Grosman,  Synthesis  of  Tschebycheff  Parameter  Symmetrical 

Filters,  Proc.  IRE,  Vol.  45,  No.  4,  pp.  454-474. (13> 

Using  Grosman's     equation  (38)  or  see  Darlington    . 
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provide  the  required  selectivity  parameter  (or  discrimination  parameters) 
and  meet  the  ripple  requirements,  if  n  is  specified  and  the  requirements 
demand  more  elements  than  is  allowed  the  computer  should  so  indicate  to 
the  operator.   It  is  estimated  that  this  information  would  be  available 
within  not  more  than  ten  minutes  from  the  start  of  computations,  and 
therefore  establish  a  doctrine  that  if  such  a  program  is  used,  the  opera- 
tor must  wait  until  this  point  in  the  program  is  passed  before  allowing 
the  computer  to  proceed  unattended,  and  by  so  doing  be  assured  that  the 
results  will  be  realizable  as  a  physical  network,  and  that  computation 
time  will  not  be  wasted. 
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5.   Conclusions 

The  results  of  the  tests  on  the  routines  programmed  to  date  indicate 
that  a  library  of  network  synthesis  routines  employing  conventional  tech- 
niques is  practicable. 

The  use  of  such  a  library  should  result  in  considerable  savings  in 
time  and  effort  for  those  people  concerned  with  network  design. 

The  accuracy  obtainable  is  more  than  sufficient  to  guarantee  excellent 
approximations  (and  undoubtedly  synthesis)  if  the  number  of  elements  to 
be  used  is  not  severely  restricted. 

The  time  required  for  solution  is  not  excessive  in  comparison  with 
the  magnitude  of  the  problem  and  the  accuracy  of  the  results.   The  employ- 
ment of  the  library  of  programs  will  require  a  working  knowledge  of  the 
computer,  unless  magnetic  tapes  are  made  up  from  the  "building  blocks" 
to  provide  specific  types  of  network  designs.   In  this  case,  detailed 
instructions  may  suffice  to  permit  persons  not  familiar  with  the  computer 
to  use  the  programs  to  good  advantage. 

The  installation  of  the  CDC  -  MK  1  Computer  will  enhance  the  value 
of  this  library  despite  the  conversion  problem.   The  Increased  speed  of 
computation  and  shorter  programs  which  will  be  possible  due  to  the  "built- 
in"  floating  point  arithmetic  will  make  such  a  library  even  more  valuable. 
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APPENDIX  I 

It  is  the  intention  of  this  appendix  to  explain  with  a  minimum  of 
theoretical  reasoning  the  overall  flow  chart  of  Ba  Hli's  method. 
The  following  inputs  are  required: 

1.  \f    (t)   -   time  sequence  of  output  function  (in  octal  or  decimal). 

2.  If  (t)(.  -   time  sequence  of  input  function  (in  octal  or  decimal), 
number  of  discontinuities. 

-  number  of  points  of  inflection, 
number  of  elements. 

6.  t      -  number  of  f.(t)  ordinates  used. 

7.  Octal  or  Decimal  input  data  indicator. 

Annex  1  shows  a  block  diagram  flow  chart  which  considers  only  inputs 
and  outputs  to  major  program  subdivisions,  and  major  logical  decisions 
required. 

We  allow  for  either  decimal  or  octal  digit  input  sequences,  in  either 
case  they  are  converted  to  binary  floating  point.   If  the  number  of  ele- 
ments (N  )  is  not  specified,  compute  AS  «  1.5  N  +  N   (see  page  17, main 
text). 

The  synthetic  division  process  is  self-explanatory,  the  quantity 
"t",  the  number  of  f.(t)  ordinates,  is  required  in  this  block. 

Annex  2  shows  a  detailed  flow  chart  for  the  determination  of  "r", 
the  relative  degree  between  numerator  and  denominator.   The  first  order 
approximation  used  to  compute  h   is  considered  to  be  sufficiently  accurate 
to  ascertain  whether  or  not  h(t)  starts  from  a  finite  value,  or  from 
zero.   The  "0"  used  in  the  decision  boxes  is  actually  slightly  greater 
than  zero,  about  .025,  to  avoid  any  erroneous  decision  that  h   is  greater 
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than  zero,  when  it  is  only  greater  due  to  the  error  introduced  by  the 
first  order  approximation.   Note  that  we  have  deviated  from  Ba  Hli's 
method  of  determining  the  behavior  of  h(t)  near  t  »  0  in  that  we  compute 
derivatives. 

We  can  now  determine,  D  ,  the  degree  of  the  denominator,  and  N  , 
the  degree  of  the  numerator  of  Z(p)  and  the  number  of  terms  of  H*(s), 
H  ,  we  must  compute  to  solve  for  Z(p). 

We  can  make  a  somewhat  arbitrary  decision  as  to  how  accurately  we 
should  compute  our  H(s),  i.e.,  what  subdivisions  we  should  use  in  our 
"moment"  computing  scheme  (see  Appendix  II),  depending  on  the  number  of 
terms  we  shall  have  in  our  denominator. 

The  next  two  appendices  give  the  program  details  for  the  solution 
of  H*(s),  and  Z'(p). 

The  solution  of  Z(p)  is  now  complete,  unless  we  desire  to  compensate 
for  dissipative  losses.   This  is  easily  accomplished  after  factoring  the 
numerator  and  denominator  polynominals.   We  preshift  our  pole  and  zero 
pattern  to  the  right  an  amount  "ex.",  unless  this  shift  causes  a  pole  or 
zero  to  move  into  the  right  half  plane  and  thereby  violate  the  conditions 
of  physical  realizability.   We  compute 
c\      -   R/2L 

2L/R        -      1/oc 

ojL/R        -    *V2^ 

2Q  /to  m        cA. 

^       »   Q/TTf 
and  thus  we  compute  <  as  a  function  of  Q  and  f ,  where  f  is  the  lowest  sig- 
nificant frequency  involved  and  is  entered  by  the  designer  as  input  data. 


fo^> 


fA(t) 
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0/D 


~V 


_. 


IIQII      &      Itfll 
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See  Appendix  II 
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APPENDIX  II 
CALCULATION  OF  H*(s)  FROM  THE  [h]    SEQUENCE 

We  are  now  required  to  form  the  power  series  expansion  of: 

Z  a  e 
i-1 

As  was  mentioned  in  the  main  body  of  this  paper,  replacing  the  Laplace' 
Stieljes  integral  by  the  summation  indicated  above,  resulted  in  an  increas- 
ing error  in  those  terms  higher  than  the  first  degree.   We  recognize  that 
this  error  results  from  the  failure  to  compute  the  exact  moment  of  |h\ 
about  the  origin. 

t"  K 

The  true  j   moment  of  an  area  defined  by  f(x)  is: 

M.  «   /   x^  f(x)  dx 
jo 

which  we  have  replaced  by  the  summation 


M.(app)  »   2 
J         i-1 


2jL  -  li  J 


2N 


A. 

1 


where  A,  *  incremental  area  or  h    and  S  ■  NK  (K«  1,2,3..,,   ) 
i  a 

For  our  program,  we  limit  K  to  4;  that  is,  we  require  the  impulse 
response  to  be  defined  over  the  range  0  to  4  time  units.   The  input  and 
output  functions  may  be  sealed  so  this  will  be  adequate  for  most  problems. 

It  might  be  well  before  proceeding  to  take  a  close  look  at  the  manner 
in  which  we  form  the  coefficients  for  the  Dirichlet  power  series  expansion. 

We  may  proceed  in  such  a  manner  as  to  form  one  term  at  a  time,  that 
is,  referring  to  Figure  II-l,  we  may  sum  first  down  the  J  »  0  column,  then 
the  j  *  1  column,  etc.   The  operation  count  for  this  procedure  is: 

N  "a"  for  the  first  term, 
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N  "a"  +  N"m"  +  "m"  for  the  second  (If  we  factor  out  the  1/2N  until 

the  sum  is  complete), 

N  "a"  +  N(j  +  l)"m"  +  (j  +  l)"m"  for  the  jth  term. 

a  total  of 

(j  +  1)  N"a"  +  N+l   2   (i  +  1)  "«" 

i-1 


(j  +  1)  JN  fa"  +  (*y.  +  ij  .V.J   +(  ^  +  l)  "m" 


where  "a"  stands  for  floating  point  addition  and  "m"  stands  for  floating 
point  multiplication. 

For  N  »  40,  j  ■  20;  this  becomes  840  "a"  +  11040  ","  and  if  we  assume 
that  floating  point  addition  requires  twice  the  time  that  floating  point 
multiplication  does,  a  equivalent  total  of  12,720  floating  point  multi- 
plications . 

If,  however,  one  proceeds  out  the  rows  and  forms  all  the  partial  sums 
associated  with  A  ,  then  A? ,  etc.,  for  i  »  1,  j  »  0,  to  j,  the  count  per 
row  is : 

"a",  2"m"  +  "a",  "a"  +  V,  "a"  +  "m" "a"  +  "m"  and  the  2N  in 

the  denominator  is  included. 

The  count  is  the  same  for  all  rows,  therefore,  the  total  is: 

N(j+l)("aW),  which  for  N  »  40,  j  -  20,  is  840  "a"  ♦  840  "m",  or 
based  on  our  previous  assumption,  a  total  of 

2,520  floating  point  operations. 
There  is,  however,  some  additional  "call-up"  and  "put-away"  address  compu- 
tation which  degrades  this  efficiency  somewhat.   Nonetheless,  it  seems 
quite  evident  that  a  saving  of  at  least  507.  can  be  effected  in  this  way. 

With  this  saving  we  can  afford  to  reduce  the  error  in  the  final  co- 
efficients by  introducing  an  interpolation  routine  which  effectively 
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doubles  or  even  quadruples  the  number  of  sample  points  used. 

The  interpolation  precedure  makes  a  straight  line  approximation  be- 
tween sample  points,  and  computes  the  value  the  intermediate  sample  points 
would  have.   Therefore,  it  is  able  to  use  a  smaller  increment  on  the 
abscissa  and  thus  approaches  a  true  integration  more  closely  than  would 
be  otherwise  possible.   The  results  of  some  runs  using  various  "degrees" 
of  interpolation  are  presented  in  Table  II-l  for  comparison  purposes. 
The  flow  chart  and  program  are  appended  hereto. 
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i-A 
2N  Al 
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(2N)     Al 
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3     2 
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2N  A3 
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2N 


u*/    \       ^  a  ^  /2N-iv    .  ^>  ,2N-i  2    .  ^>  /2N-lv3    A  ^  /2N-lv-i    . 

H*(s)   -|Ai+2(-2F)   Ai+f  (1F>     Ai+|("2F)     Ai+    '••   +f  <~ >     Ai 

The  factorials  which  appear  in  H(s)  are  yet  to  be  supplied,  and  are 
brought  in  with  the  alternating  sign  according  to  whether  j  is  odd  or  even. 


Formation  of  R*(s) 
Figure  II-l 
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DIRICHLET  POWER  SERIES  EXPANSION  FOR  PSEUDO-SYSTEM  FUNCTION 
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17 

2010 

3000 

0305 

If  TS  1  up,  0305  (l) 

0305 

30 

0556 

2100 

05^5 

No,  Set  t  to  1 

030^ 

31* 

3000 

2100 

0311 

and  skip  to  03 11 

0305 

17 

2020 

3000 

0310 

If  TS  2  up,  0310 

0306 

30 

05  *u 

2100 

05^5 

No,  Set  t  to  2 

0307 

3U 

3000 

2100 

0311 

and  skip  to  0311 

0310 

30 

05^2 

2100 

0545 

Set  t  to  3 

0311 

35 

0536 

05^5 

0555 

N  exp  +  t  to  N'  exp 

0312 

32 

0510 

0516 

031^ 

Set  h*  pickup 

0313 

31 

0533 

0574 

0531 

"-1"  to  SA 

031U 

31 

0120 

0121 

0523 

ht   to  TF 

0315 

3k 

05^5 

0532 

0321 

If  t  *  1,  0321 

0316 

31 

0523 

0564 

0530 

TF  to  TA 

0317 

30 

0553 

2100 

05V7 

1  to  H 

0320 

3* 

3000 

2100 

0413 

Skip  to  OU13 

0321 

3k 

05^5 

0534 

0350 

If  t  *  2,  0550 

0322 

33 

0572 

0533 

0326 

No,  if  SA  *  "0",  0326 

0323 

36 

0525 

0533 

0521 

No,  TF/2  to  A 

032^ 

30 

0564 

2100 

0562 

0325 

3h 

3000 

2100 

0333 

Skip  to  0333 

0326 

31 

0523 

0564 

2000 

Load  TF 

0327 

31 

052U 

0565 

2002 

Load  TE 

0330 

3k 

3000 

2100 

1516 

TF  -  TE 

0331 

36 

2000 

053^ 

2000 

(TF  -  TE)/U 

0332 

31 

2000 

2001 

0521 

(TF  -  TE)/U  to  A 

0335 

31 

0523 

05  6k 

2000 

Load  TF 

03  3U 

31 

0521 

0562 

2002 

Load  6 

0335 

3k 

5000 

2100 

1516 

TF  -  A 

0336 

36 

2000 

0532 

2000 

(TF  -A)/2 

0337 

31 

2000 

2001 

0530 

to  TA 

03^0 

51 

0523 

05  6k 

2000 

Load  TF 

O3U1 

31 

0521 

0562 

2002 

Load  & 

03^2 

3^ 

3000 

2100 

15^0 

TF  4-  A 

03^3 

36 

2000 

0532 

2000 

(TF  +  A)/2 

03^ 

31 

2000 

2001 

0527 

to  TB 

03^5 

31 

0523 

0564 

052U 

TF  to  TE 

05^6 

30 

05  3k 

2100 

05^7 

2  to  H 

03^7 

3k 

3000 

2100 

OU13 

Skip  to  04l3 

0350 

33 

0572 

0533 

0356 

Is  "SA"  *  "0" 

0351 

36 

0523 

053^ 

0521 

No,  TF/U 

0352 

30 

0564 

2100 

0562 

to  A 

0353 

56 

0523 

0532 

0522 

TF/2 

035^ 

50 

056U 

2100 

0563 

to<^ 

0355 

3k 

5000 

2100 

0565 

Skip  to  O365 

0356 

31 

0523 

056^ 

2000 

Load  TF 

0357 

51 

052U 

O565 

2002 

Load  TE 

(1 )  The  test  switch  function  would  be  performed  by  the  interpolation 
selection  in  the  overall  synthesis  program. 
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DIRICHLET  POWER  SERIES  EXPANSION  FOR  PSEUDO-SYSTEM  FUNCTION 


ADD       INST 


MJ 


M* 


M3 


REMARKS 


O36O 
O36I 
O362 
0363 
03  6k 
0365 
O366 

0367 
0370 
0371 
0372 
0373 
037^ 
0375 
0376 
0377 

01*00 
01*01 

01*02 
0U03 
0t*0U 
0U05 
Ol*06 


OUll 
0U12 
OU13 
OUlU 

01*15 
01*16 
0U17 
01*20 
01*21 


0U2U 
01*25 


0430 
01*31 
01*32 
01*33 

01*3** 
0^35 
0^36 

0U37 


3* 
36 
31 
36 
31 
31 
51 
3^ 
31 
3U 
36 
30 
31 
3* 
56 
30 
31 
31 
3^ 
36 
50 
31 
3^ 


01*07   36 

01*10  30 


31 
30 
31 
31 

3^ 
31 
31 
3U 
31 


0U22   32 
0U23   33 


35 
3k 


0U26  3k 
0l*27   31 


31 
3k 
31 
31 
31 
3k 
30 
32 


3000 

2000 

2000 

2000 

2000 

0523 

0522 

3000 

0521 

3000 

2000 

2001 

0522 

3000 

2000 

2001 

0523 

0521 

3000 

2000 

2001 

0522 

3000 

2000 

2001 

0523 

05^3 

0531 

05  3^ 

3000 

2000 

0535 

3000 

2000 

0460 

0520 

2100 

05  M* 

3000 

0530 

0520 

3000 

2000 

0600 

0530 

3000 

0U33 

2006 


2100 

053^ 
2001 
0532 
2001 
0564 
0563 
2100 
0562 
2100 

05  3^ 
2100 

0563 
2100 

053^ 
2100 
05  6k 
0562 
2100 

053^ 
2100 

0563 
2100 
05lk 
2100 
0564 
2100 
0572 

0575 
2100 
2001 
0576 
2100 
2001 
0516 
2100 
2100 
2100 
2100 
0571 
0561 
2100 
2001 
0601 
0571 
2100 
1505 
1677 


1516 
2000 
0522 
2000 
0521 
2000 
2002 
1516 
2002 
1516 
0530 
0571 
2002 
15U0 
0527 
0570 
2000 
2002 
15U0 
0526 
0567 
2002 
151*0 
0525 
O566 
052U 
05  V7 
2000 
2002 
151*0 

0551 
2001* 

1507 
0520 
01*33 
01*57 
051*1* 
01*27 
01*33 
2000 
2001* 
1500 
0530 
2000 
2002 
151*0 
2006 
01*1*2 


TF  -  TE 
(TF  -  TE)/1* 

to  * 
(TF  -  TE)/8 

to  A 
Load  TF 
Load  <*. 
TF  -  <* 
Load  A 
TF  -  «K  -  A 
(TF  -  *  -  A)A 

to  TA 
Load  A 
TF-  A 
(TF  -  &)/k 

to  TB 
Load  TF 
Load  A 
TF  +  A 
(TF  +  A  )/U 

to  TC 
Load  <* 
TF  +A  +  * 
(TF  +A+0O/I* 

to  TD 
TF  to  TE 
Set  H  to  1* 
Load  SA 
Load  "2" 
SA  +  2 

to  SA 
Load  N1 
Divide 
Result  to  M 
Reset  putavay 
If  M  *  0,  01*57 
No,  clear  J  tally 
If  J^O,  skip  to  01*27 
No,  skip  to  01*33 
Load  TA 
Load  M 
TA  X  M 

to  TA 
Load  summation 
Load  TA 
Accumulate 
Build  putavay 
Build  putavay 
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DIRICHLET  POWER  SERIES  EXPANSION  FOR  PSEUDO-SYSTEM  FUNCTION 


ADD 

INST 

M1 

M* 

M3 

0440 

35 

0442 

0532 

2006 

044l 

32 

2006 

1677 

0443 

0442 

30 

2000 

2100 

06l6 

0443 

30 

2001 

2100 

0617 

0444 

35 

0544 

0532 

0544 

OUU5 

35 

0^33 

0512 

0453 

0446 

34 

05  40 

0544 

0425 

0447 

36 

05U7 

0532 

0547 

0450 

3^ 

05U7 

2100 

0453 

OU51 

35 

031^ 

0512 

0314 

OU52 

3^ 

3000 

2100 

0314 

OU53 

31 

0527 

0570 

0530 

0454 

31 

0526 

0567 

0527 

OU55 

31 

0525 

0566 

0526 

OU56 

3^ 

3000 

2100 

0413 

0457 

32 

05  14 

0516 

0464 

0460 

31 

06OO 

0601 

1310 

046l 

34 

3000 

2100 

1011 

0462 

30 

0532 

2100 

0544 

0U63 

31 

0532 

0573 

0517 

o464 

31 

0620 

0621 

2000 

0465 

32 

0544 

0532 

2106 

0466 

3^ 

2006 

2100 

0470 

0^67 

3^ 

3000 

2100 

0471 

0470 

32 

044l 

0537 

2001 

0471 

31 

0517 

0560 

2004 

0U72 

3^ 

3000 

2100 

1507 

0473 

31 

2000 

2001 

1310 

0474 

3^ 

3000 

2100 

1011 

0475 

35 

0544 

0532 

0544 

0U76 

31 

0551 

0544 

2000 

0477 

31 

0517 

0560 

2004 

0500 

3^ 

3000 

2100 

1500 

0501 

31 

2000 

2001 

0517 

0502 

35 

0464 

0512 

0464 

0503 

3^ 

0540 

0544 

0464 

0504 

22 

0000 

0000 

0000 

0505  • 

■507 

Available 

0510 

00 

0000 

0001 

0000 

0511 

35 

2000 

2005 

2000 

0512 

00 

0002 

0002 

0000 

0513 

34 

3000 

2100 

0302 

0514 

00 

0602 

O603 

0000 

0515 

00 

0000 

0000 

0001 

0516 

00 

7777 

7777 

0000 

0517  - 

•  0556 

Used  for  tempoary  stc 

0537 

02 

0000 

0000 

0000 

0540 

00 

0000 

0000 

0000 

REMARKS 


Build  put away 
Build  put away 
Putaway 
Putaway 

Add  1  to  J  tally 
Modify  summation  pickup 
If  J  *  J,  repeat  to  0425 
No,  reduce  H  by  1 
If  H  •*-  0,  skip  to  0453 
No,  modify  h±   pickup 
Repeat  loop 
Shift  TB  to  TA 
TC  to  TB 
TD  to  TC 
and  repeat  loop 
Reset  hi  pickup 
First  terra  pickup 
Convert  and  print  HQ  (2) 
Set  J  to  1 
Load  F! 

Pickup  Hj  term 
Extract  last  bit  of  J 
If  J  odd,  skip  to  0470 
No,  skip  to  0471 
Change  sign  to  minus 
Load  F! 
Divide  by  F! 
Load  H 

for  conversion  and  print  out1 
Increment  J  by  1 
Load  "J" 
Load  F! 

Multiply  for  (F  +  l)« 
Putaway 

Modify  H±   pickup 
Have  J  terms  been  divided  by 

appropiate  factorial;  no,  0464 
Yes,  halt 

h^  pickup  reset 
Putaway  clear  routine 
Modifier  for  0502 
Putaway  clear  routine 
Putaway  reset 
Putaway  clear  routine 
M1  and  M2  extractor 

Sign  extractor 
J 
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DIRICHLET  POWER  SERIES  EXPANSION  FOR  PSEUDO-SYSTEM  FUNCTION 


ADD 

INST 

M1 

M2 

M3 

REMARKS 

05^1 

00 

0000 

0000 

0002 

2  -  for  H 

05^2 

00 

0000 

0000 

0003 

3  -  for  H 

05^3 

00 

0000 

0000 

000U 

k   -  for  H 

Ojkk 

00 

0000 

0000 

0000 

J  tally 

05^5 

00 

0000 

0000 

0000 

t  storage 

05^6 

00 

0000 

0000 

0001 

1 

05^7 

00 

0000 

0000 

0000 

H  storage 

0550 

26 

2100 

2100 

0600 

Put away  clear  routine 

0551 

00 

0000 

0000 

OO^U 

Constant  for  OU76 

0552 

3h 

200U 

2000 

2000 

Putaway  clear  routine 

0553 

00 

0000 

0000 

0000 

Available 

055^ 

00 

2100 

2100 

1000 

Testword  for  clear 

0555  ■ 

-  0577 

Used 

for  tempo 

rary  stora 

ge 

0600  -  1000  Used  for  storage  of  results  (may  be  reduced  to  provide 
a  minimum  of  2 (J  +  l)  cells. 


(2)  The  conversion  and  print  routine  is  linked  in  this  program  for 
test  purposes  only  and  should  be  replaced  with  a  boot -strapping 
routine  to  call  in  next  portion  of  the  overall  program. 
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FLOW  CHART  OF  DIRICHLET  POWER  SERIES  EXPANSION 


Set  K-^-4 


Kr- 


JL 


Clear  Tally 


Initialize  B  test  Command 


Add  1  to  tally 
B  -  2->B 


No 


Yes 


K-l  — >K 


JL 


K  -  1-»K' 

Clear  Putaways 
Set  h  pickup 


N    +  t-^.  N' 
exp         exp 

"-1" >   SA 


-^ 


J 


©- 


H-l  -»  H 


Yes 


mod  h     pickup 


o 


_  J 


SA  +   "2" 
SA  +     N' 
Reset  PA 


Yes    /M>K'\  No 


o-*j 


No 


0425 


Yes 


IA  x  M-*TA 


<]/    0433 


Set  H  pickup 
Pickup  H 
Convert  &  Print 
Set  j  -  1 


^F! 


Pickup  Ht-»2000 
Extract  last  bit 

°f  i 


TB 
TC 
TD 


TA 
TB 
TC 


Yes 


ex(-)  into  2001 


Div  by  Fi 
Conv  &  Prt 
add  1  to  j 
j_^  „.,„ 

"J»x  FJ-^  F! 

mod  H  pickup 
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APPENDIX  III 
MATRIX  INVERSION 

As  was  noted  on  pages  14  and  15  of  the  main  text  of  this  paper,  the 
final  solution  of  Ba  Hli's  method  led  to  a  system  of  simultaneous  linear 
equations.   While  it  is  certainly  possible  to  program  the  solution  of  a 
system  of  equations  by  the  straight  forward  elimination  technique,  this 
method  is  awkward  when  many  equations  are  involved  (and  for  precise 
synthesis  of  transcendental  functions,  there  must  be  many).   The  computer 
must  do  all  the  address  computation  based  on  N,  the  number  of  poles  (or 
the  degree  of  the  denominator)  and  r,  the  relative  degree  between  numera- 
tor and  denominator.   This  method  for  n  equations  in  n  unknowns,  neglect- 

2 
ing  the  address  computations,  would  require  n  divisions  to  eliminate  the 

2 
first  unknown,  (n-1)   to  eliminate  the  second,  etc.   In  all,  to  reduce 

the  system  to  the  desired  "one  equation  in  one  unknown"  (we  are  dealing 

with  square  matrices),  requires: 


£   .2   n(n+l)(2n+l)       _    . 
/       l   ■   *    '} *•  operations* 


i-1 

If,  and  only  if,  we  have  saved  the  intermediate  equations,  we  may 

complete  the  solution  in 

2 
n  +n       ... 
— - —  operations 

or  for  the  complete  solution, 

3    o   o 
n     2    2n 

-r-  +  n  +  —t     operations 

Such  a  procedure  would  require  an  inordinate  amount  of  storage 


*The  term  "operation"  is  customarily  taken  to  mean  multiplications 
or  divisions.  Addition,  subtractions  and  address  modifications  are  not 
included  in  operation  counts. 
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space  and  worse,  intricate  address  computations. 

For  simplicity  of  programming  and  best  relative  speed,  we  would 
suggest  that  the  denominator  matrix,  i.e.,  that  group  of  equations  which 

are  equated  to  zero,  and  from  which  we  determine  the  q's,  be  solved  by 

1  2 
the  elimination  form  of  the  inverse,  '   and  solve  the  remaining  equa- 
tions for  pfl,  p. p  ,  which  incident  tally  will  always  be  explicit  in 

terms  of  the  q's,  by  simple  algebra. 

We  will  describe  this  highly  efficient  method  of  computing  an  in- 
verse matrix  and  demonstrate  its  simplicity  and  economy  of  operation 
by  means  of  a  simple  example.   Briefly,  by  way  of  introduction,  we  have 
the  situation: 

[H]    [Q]    »    [h]  ,    if  we  multiply  by  the   inverse    [h]"1 

[H]        [H]     [Q]    «    [H]"1    [h]      we  have 

[Q]  -    [k]"1    [hi,     the  desired  result. 

Our  immediate  problem  is  to  find  a  rapid  means  of  computing  the 
inverse  of  H.   The  general  method  of  solution  can  be  summarized: 

(1)  [H  l]  This  operation  can  be  carried  out  by  our 

(2)  [H   H  H  '  ij    suggested  technique  operating  solely  on 

(3)  [I  H"1]  the  unit  matrix. 

We  perform  this  transformation  on  I  by  bringing  in  H  one  column  at 
at  time,  multiplying  by  those  elements  which  have  been  transformed  by  the 
introduction  of  the  columns  preceding. 

George  B.  Dantzig,  Computational  Algorithm  of  the  Revised  Simplex 
Method,  USAF  Project  Rand  Memorandum,  ASTIA  Document  No.  AD  114135,  26 
Oct.,  1953. 

2 
H.  M.  Markowitz,  The  Elimination  Form  of  the  Inverse  and  Its 

Application  to  Linear  Programming,  Management  Science,  Vol.  3,  No.  3, 

April,  1957. 
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The  actual  transformation  may  be  summarized 


(1)  H 

(2)  I 


j ,k  a  1,2  ...  .n 


<3>  Vk-1^ 

<4>     X]k  ■  V  *  Pkr  «   -  r> 

<5>     Xjk  *  V  "  Pkj    lk  (^r) 


1 , 2  .  .  .  .n 


which  states  (for  those  not  versed  in  matrix  notation),  given  a  n  x  n 
matrix  H  (1),  construct  the  corresponding  n  x  n  unit  matrix  I  (2). 
Multiply  the  k —  column  of  H  by  the  unit  matrix  (or  the  partly  trans- 
formed matrix)  to  produce  the  pivotal  column  P   ,  (3).   Divide  the  ele- 

t  H  #*Vi 

ments  of  the  r —  row  of  the  unit  or  matrix  under  transformation  by  r — 
element  of  the  pivotal  column,  (4).   Correct  all  other  rows  by  subtracting 
from  the  element  to  be  corrected  the  product  of  pivotal  element  in  that 
row  and  the  corresponding  element  in  the  i —  row  computed  in  Step  (4), 
(5). 

A  simple  example  will  make  the  procedure  clear,  say  H  is 

4   2   3 

10  2 

2   12 
We  start  by  writing  down  the  unit  matrix  of  corresponding  order,  viz., 


(2) 


K  i 
r 

(3) 


10  0 
0   1   0 

1,  R  -  1 
"10  0  4 
0   10   1 
0  0   12 


Note:   This  is  akin  to  merely  augment- 


ing the  unit  matrix  by  H   for  R  =  1. 


Ill 


K  =  1 

r 

,  R  - 

1 

(4) 

"1/4 

0 

0 

4 

0 

1 

0 

1 

0 

0 

1 

2 

(5) 

"1/4 

0 

0 

• 

-1/4 

1 

0 

• 

-1/2 

0 

1 

• 

lL  -  xik  *  \- 

j   -  r  -   1 


1,2. . .  .n 


l'n  '  ln  -  pu  xik  (J  < '  -  D<k  -  i) 

Discard  P.  .   Note:   Columns  to  the  right 

j 

of  I  ,   are  not  operated  on,  herein  is  the 

^  r 
great  economy  of  this  method. 


r 
(3) 


(4) 


(5) 


2,  R  • 

"1/4 
-1/4 
-1/2 
"1/4 

1/2 
1-1/2 
"  0 

1/2 
-1/2 


2 

0  0 

1  0 
0  1 
0  0 

-2  0 

0  1 

1  0 
-2  0 

0  1 


1/2 
■1/2 
0 

1/2 
■1/2 
0 


V  •  H2 


J'   =  J„   +  P„   (second  element  in  column  P) 


Note:   Again  since  K  =  2<3,  column  3  was 


unaffected;  also  note  since  P, 


was  unaffected. 


0,  I 


3k 


K  »  3 

,  our  final 

trans  formation 

is : 

(3) 

0   1   0 

2  ~ 

(4) 

~0 

1 

0 

2~ 

(5) 

"2 

1 

-4 

-1/2  -2   0 

-5/2 

1/2 

-2 

0 

-5/2 

-4/2 

-2 

5 

-1/2   0   1 

1/2. 

:1 

0 

2 

l/2_ 

_"1 

0 

2 

and  therefore 

~2 

1   -4 

H"1    . 

-2 

-2   5 

-1 

0 

2 
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Check 


1 


M  DO"  -DO 

4   2   3 

1  0   2 

2  1   2 


2 

1 

-4 

-2 

-2 

5 

-1 

0 

2 

The  total  number  of  operations  (neglecting  all  incidental  nulls)  is: 


Multiplications  to  bring  in  each  new  row 

Multiplications  and  divisions  to  accomplish 
each  transformation 

To  solve  the  equations  once  the  inverse  has 
been  computed 

Total  product  operations  to  complete  solution 


n(n+l)(n+2)/3 

(n+l)(n)(n-l)/6 
2 


3     2 

n  +  4n  +  n 


Reference  to  Table  III - 1  shows  the  efficiency  of  the  elimination 
form  of  the  inverse.   It  is  important  to  note  that  by  minimizing  the 
number  of  operations,  we  also  minimize  the  propagation  of  round-off  error. 

Von  Neumann  has  determined  that  an  approximate  inverse  will  be  found 


if 


s/4 


n  <  .  15 /B    ,  where  G>   is  the  base  of  the  number  system  used  in  the 
computation,  and  S  is  the  number  of  places.   The  program  for  the  elimina- 
tion form  of  the  inverse  (appended)  used  a  packed  floating  point  where 
the  word  structure  was: 

XS  MMMM  MMMM  MQEE  Legend:   X  -  not  used 

S  ■  sign  of  magnitude 

M  ■  magnitude  (octal  digit) 

Q  «  2  binary  bits  of  magnitude 

&  sign  of  exponent 
E  a   magnitude  of  exponent 

J.  Von  Neumann  and  H.  H.  Goldstine,  Numerical  Inverting  of  Matrices 
of  High  Order,  American  Mathematical  Society  Bulletin,  Vol.  53,  pp.  1096. 
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This  realized  an  S  of  29,  and  therefore  theoretically  the  program  should 
successfully  compute  the  inverse  for  all  matrices  of  n<23. 

The  maximum  round -off  error  to  be  expected  is  (after  Von  Neumann): 

which  f or  Q «  2,  S  ■  29,  and  n  «  25  reduces  to 

lfi  1  2  S 
€  ■   228   "  -000067 

Thus  we  can  expect  an  accuracy  of,  at  worst,  four  significant  figures  for 

a  25  X  25  matrix. 

If  we  decide  to  use  a  double  precision  floating  point,  S  (for  the 

CRC  102A)  would  increase  to  66,  by  using  the  following  configuration: 

XS  .  MMMM  MMMM  MMMM      XS   MMMM  MMMM  MMEE 

where  S   is  sign  of  magnitude 

S   is  sign  of  exponent 

and  the  round -off  error  now  becomes  (for  n  »  25) 

1.8125  ^    1ft-9 

*    '  3.7  X   10U>    K<    10 

This  requires  a  considerably  greater  amount  of  computation  time 

since  a  floating  point  multiplication  is  now  three  "old"  floating  point 

multiplications  and  two  floating  point  additions,  however  even  for  n  »  60, 

-9 
6  <  10   .   For  comparative  purposes,  note  Table  III-2.   The  choice  of 

either  packed  or  double  precision  is  really  dependent  on  the  order  of 

the  matrix  expected  and  the  requirements  of  the  problem,  it  seems  to  this 

writer  that  both  procedures  should  be  available. 

It  roust  be  mentioned  in  passing  that  the  method  of  Hotelling  and 

Bodewig   for  iterating  to  find  an  improved  inverse  from  an  approximate 

A.  S.  Householder,  Principles  of  Numerical  Analysis,  McGraw-Hill, 
New  York,  1953,  pp.  80  ff. 
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inverse  and  the  original  matrix,  is  readily  available.   However,  although 

the  programming  of  this  scheme  would  be  extremely  simple,  inasmuch  as 

3 
the  basic  matrix  multiplication  is  already  present,  two  n  multiplica- 
tions are  required  for  each  successive  improvement,  and  furthermore, 
another  n  x  n  matrix  must  be  stored.   These  requirements  with  regard  to 
storage  and  computation  time  make  the  double  precision  floating  point 

arithmetic  routine  more  attractive,  since  the  time  to  solve  an  n  x  n  and 

5n3   4n2  -  n 
perform  one  iteration  would  be  proportional  to  r X  200  for 

the  packed  floating  point,  whereas  with  equal  storage  space,  the  double 

n3   4n2  +  n 
precision  would  only  be  proportional  to  - X  425. 

The  decision  as  to  whether  or  not  the  simple  packed  floating  point 
routine  or  the  double  precision  was  to  be  used  would  be  made  in  the  initial 
stages  of  the  problem;  i.e.,  when  "n  +  r"  was  determined,  assuming  the 
number  of  elements  was  not  specified  (see  Appendix  I). 

In  order  to  obtain  a  practical  estimate  of  the  accuracy  of  the 
matrix  inversion  program  and  the  time  required,  the  following  tests  were 
made : 

A  symmetric  10  x  10  matrix  was  formed  from  a  column  of  random  numbers, 
the  sign  and  position  of  the  decimal  point  were  also  established  using 
random  numbers.   This  matrix  was  then  inverted  using  the  program.   The 
matrix  product  of  the  original  matrix  and  the  inversion  was  then  com- 
pared with  a  unit  matrix.   The  rms .  error  of  each  element  was  found  to 
be  431  x  10    ,  the  maximum  error  was  2430  x  10 

A  similar  test  for  6x6  matrix  resulted  in  an  rms.  error  of  14  x 

-9  -9 

10   and  a  maximum  error  of  43  x  10 

3 
The  running  time  was  determined  to  be  approximately  .064  n  minutes, 

(64  minutes  for  a  10  x  10,  14  minutes  for  a  6  x  6). 
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The  flow  chart  and  program  for  the  elimination  form  of  the  inverse 
using  simple  packed  floating  point  airthmetic  is  appended  hereto. 
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TABLE  III-l 

COMPARISON  OF  VARIOUS  METHODS 

Method  No.  of  Operations 

1  3 

Stiefel-Hestenes  2n  +  .... 

1  3    2 

Strict  Inversion  2n  +  n 


Orthogonalization 


1  2n3  ♦  „2/2 


2  3     2 

Elimination  Form  1/2  (n  +  2n  +  n) 

3     2 
Back-Substitution  1/6  (2n  +  5n 


♦Assumes  product  operation  requires  70  word-times. 


Word -Time 

!S  (n> 

-20)* 

112 

X 

io4 

115 

X 

IO4 

113 

X 

io4 

31 

X 

io4 

21 

X 

io4 

TABLE    II I -2 
COMPUTATION  TIME   IN  WORD-TIMES 


Operation 

Normal  Float 

Packed  Float 

Double  Precision  Float 

Multiply 

80 

200 

425 

Divide* 

112 

232 

684 

♦Assuming  overflow  occurs  1/4  of  the  time. 


Householder, 

2x,       x,       (24) 
Von  Neumann , 
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MATRIX  INVERSION 

PROGRAM 

(CRC  -  102A) 

ADD 

INST 

M1 

M2 

M3 

REMARKS 

0100 

30 

0076 

2100 

0371 

Store  N 

0101 

30 

0077 

2100 

0366 

Store  (n  +  r) 

0102 

30 

0371 

0561 

0355 

Shift  N  to  M1  position 

0103 

05 

3000 

3000 

0370 

Clear  cells 

010U 

5»* 

3000 

2100 

2000 

0500  -  1777 

0105 

35 

0377 

2100 

0106 

Initialize  0106 

0106 

30 

0077 

2100 

1777 

Shift  the 

0107 

35 

0106 

0376 

0106 

contents  of  Channel  0 

0110 

3U 

0375 

0106 

0106 

to  Channel  17 

0111 

35 

2100 

2100 

0373 

Set  RC  (row  count)  to  0 

0112 

35 

0371 

0335 

2006 

N  ♦  1  modifier 

0113 

32 

0337 

0332 

OllU 

Initialize  OllU 

OllU 

30 

0367 

2100 

051U 

Build 

0115 

35 

OllU 

2006 

OllU 

N  x  N 

0116 

35 

0373 

0335 

0373 

unit 

0117 

3* 

0371 

0373 

OllU 

matrix 

0120 

30 

0335 

2100 

0336 

Set  CM  (column  count )  to  1 

0121 

35 

03  5U 

O366 

O565 

1701  +  (n  ♦  r)  =  S 

0122 

3U 

3000 

2100 

012U 

Skip  0123  initially 

0123 

35 

0336 

0335 

0336 

Increment  CM  by  1 

012U 

32 

2100 

0332 

0152 

Set  putaway  to  0000 

0125 

35 

2100 

2100 

0373 

Set  RC  to  0 

0126 

32 

0362 

0365 

0135 

Reset  pickup  for  ajj 

0127 

36 

0365 

0356 

0327 

S  -  CM  *  RS(  row  selector) 

0130 

26 

2100 

2100 

0323 

Clear  CC  (column  tally)  and  TS 

0131 

5U 

0327 

0326 

0133 

If  RS  greater  1677,  skip  to  133 

0132 

3U 

3000 

2100 

01 U5 

No,  skip  tp  01U5 

0135 

50 

0327 

O361 

2006 

Shift  RS  to  M1  position 

013U 

32 

2006 

0363 

01U0 

Intalize  OlUO  to  RS 

0135 

30 

0510 

2100 

2000 

Pickup  an  to  "A" 

If  A  4   0,  skip  to  OlUO 

0136 

5U 

2000 

2100 

01U0 

0137 

3U 

3000 

2100 

01U5 

If  A  e  0,  skip  to  OIU5 

OlUO 

30 

1704 

2100 

2001 

Load  column  element 

OlUl 

3U 

5UOO 

2100 

0U00 

Transfer  to  multiply 

0lU2 

30 

O36U 

2100 

2001 

Load  TS 

O1U5 

51* 

5000 

2100 

0U00 

Transfer  to  add 

OlUU 

30 

2000 

2100 

O36U 

Store  sum  in  TS 

01U5 

35 

0323 

0335 

0323 

Increment  CC  by  1 

01  U6 

35 

0135 

0360 

0135 

Increment  a^ «  pickup  by  1 
Increment  Rs  by  1 

01U7 

35 

0327 

0335 

0327 

0150 

3U 

0371 

0323 

0131 

If  N  greater  than  CC,  repeat 

0151 

36 

0135 

O36O 

0135 

No,  compensate  pickup 

0152 

30 

03  6U 

2100 

0000 

Putaway  row-column  product 

0153 

35 

0373 

0335 

0373 

Add  1  to  RC 

015U 

3* 

0371 

0373 

0156 

If  N  greater  than  RC,  skip  0155 

0155 

3* 

3000 

2100 

Ol6l 

No,  skip  to  0155,  thence  to  0l6: 

0156 

35 

0135 

O36O 

0135 

Add  1  to  a^j  pickup 
Increment  putaway  address 

0157 

35 

0152 

0335 

0152 

0160 

3* 

3000 

2100 

0127 

Repeat  loop 
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MATRIX  INVERSION 

PROGRAM 

(CRC  -  102A) 

ADD 

INST 

M1 

M2 

M3 

REMARKS 

Ol6l 

36 

0556 

0335 

0357 

CM  -  1  x  J 

0162 

26 

0357 

0371 

2006 

J  x  N 

OI65 

35 

2006 

0337 

2007 

0500  +  (J  x  N) 

Ol64 

35 

2100 

2100 

0525 

Clear  CC 

OI65 

30 

2007 

0550 

0356 

Initialize  Ir, 

0166 

32 

0356 

0563 

0172 

In  0172 

OI67 

30 

0357 

0330 

2006 

Initialize  Pk 

0170 

52 

2006 

0365 

0173 

in  0175    r 

0171 

32 

2007 

0332 

0175 

Initialize  I£  in  0175 

0172 

30 

0510 

2100 

2000 

Load  Ir 

0173 

30 

0002 

2100 

2001 

Loaci  Pk 

0174 

3U 

3600 

2100 

0400 

Ir  T  Pk  -  Ir 

0175 

30 

2000 

2100 

0510 

Put away  result 

0176 

35 

0323 

0355 

0323 

Increment  CC  by  1 

0177 

31* 

0536 

0325 

0201 

If  CM  greater  than  CC,  0201 

0200 

34 

3000 

2100 

0204 

If  CM  =  CC,  0204 

0201 

35 

0172 

0360 

0172 

Increment  I_.  to  1^  , 
A      i +1 
Increment  putaway 

0202 

35 

0175 

0335 

0175 

0203 

3* 

3000 

2100 

0172 

Continue  dividing  rov  »r* 

0204 

35 

2100 

2100 

0325 

Set  M  to  0 

0205 

3* 

0336 

0325 

0207 

If  CM  greater  than  M,  0207 

0206 

34 

3000 

2100 

0242 

If  CM  =  M,  0242 

0207 

55 

0337 

0325 

0322 

Build  I' 
Build  Iij 

0210 

30 

0322 

0350 

0324 

0211 

55 

0571 

0335 

0327 

N  +  1  to  TS 

0212 

50 

0555 

2100 

0375 

Initialize  RC  to  1 

0213 

32 

0356 

0363 

0222 

Initialize  Iir  pickup 

0214 

32 

2100 

0365 

0223 

Initialize  Pj,  pickup 
Initialize  I*[,  pickup 
Initialize  Ijj  pickup 

0215 

32 

0324 

0365 

0226 

0216 

32 

0322 

0332 

0230 

0217 

34 

0373 

0336 

0222 

If  RC  4     CM  ' 

0220 

34 

0356 

0375 

0222 

skip  to  0222 

0221 

3^ 

3000 

2100 

0231 

If  RC  »  CM,  skip  to  0251 

0222 

30 

0510 

2100 

2000 

Load  Ilr 

0225 

30 

0002 

2100 

2001 

Load  Pk 

0224 

3U 

3400 

2100 

0400 

Multiply" 

0225 

30 

2000 

2100 

2001 

Shift  for  subtraction 

0226 

50 

0510 

2100 

2000 

Load  Iij 

I'i«  -  IX1  -  Iir  x  Pk 
Putaway  J 

0227 

3U 

5200 

2100 

0400 

0230 

50 

2000 

2100 

0510 

0231 

35 

0373 

0335 

0373 

Increment  RC  by  1 

0232 

5^ 

0327 

0373 

0236 

If  N  ♦  1  greater  RC,  0236 

0253 

35 

0356 

O360 

0356 

Increment  Ilr  to  Igr 

02  34 

35 

0325 

0335 

0325 

Increment  M  by  1 

0235 

3^ 

3000 

2100 

0205 

Repeat  loop 

0236 

55 

0223 

O36O 

0223 

Increment  rV  to  Pk  . 
Iij  +  N  =  I2r     r+i 
Increment  putaway 

0257 

35 

0226 

0355 

0226 

02U0 

35 

0230 

0371 

0230 

0241 

3h 

3000 

2100 

0217 

Repeat  loop 

0242 

3k 

0371 

0336 

0123 

Next  iteration  (  «.  loop) 
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MATRIX  INVERSION  PROGRAM   (CRC  -  102A) 
ADD   INST    M1      M2      M3  REMARKS 


0243 

35 

2100 

2100 

0336 

02  44 

35 

0347 

2100 

0152 

O2U5 

30 

0367 

2100 

0000 

0246 

32 

0376 

0332 

0152 

0247 

35 

0354 

2100 

0155 

0250 

3* 

3000 

2100 

0125 

0251 

35 

031*6 

2100 

0152 

0252 

35 

0353 

2100 

0155 

0253 

35 

2100 

2100 

0352 

0254 

52 

0351 

0332 

0276 

0255 

34 

3000 

2100 

0260 

0256 

35 

0552 

0435 

0352 

0257 

55 

0276 

0455 

0276 

0260 

32 

0551 

0563 

0267 

0261 

35 

«^_100 

2100 

0364 

0262 

30 

0352 

036l 

2006 

O263 

32 

2006 

0363 

0270 

0264 

54 

3000 

2100 

0267 

0265 

35 

0267 

0360 

0267 

0266 

56 

0270 

O360 

0270 

0267 

30 

1701 

2100 

2000 

0270 

50 

0000 

2100 

2001 

0271 

34 

3400 

2100 

ouoo 

0272 

30 

0364 

2100 

2001 

0273 

34 

3000 

2100 

0400 

0274 

30 

2000 

2100 

0364 

0275 

3* 

0270 

0350 

0265 

0276 

30 

036U 

2100 

1001 

0277 

54 

0366 

0352 

0256 

0300 

35 

2100 

2100 

0345 

0301- 

0321 

Available  for 

boot-81 

Denominator  in  Chann* 

Numerator  in  Channel 

0322 

00 

0000 

0000 

0502 

0323 

00 

0000 

0000 

0003 

032U 

00 

0502 

0000 

0000 

0325 

00 

0000 

0000 

0003 

0326 

00 

0000 

0000 

1677 

0327 

00 

0000 

0000 

1705 

0330 

00 

0000 

0000 

0030 

0331 

35 

2000 

2005 

2000 

0332 

00 

0000 

0000 

3777 

0333 

34 

3000 

2100 

0105 

0334 

00 

0000 

0000 

1701 

0335 

00 

0000 

0000 

0001 

0336 

00 

0000 

0000 

0000 

0337 

00 

0000 

0000 

0500 

Set  CM  =  0 

Reset  putavay 

Putavay  "l"  in  0000 

Dummy  (see  0244) 

Mod  0155  ("out  comd") 

Transfer  for  N+l8  column 

Reset  putavay 

Restore  0155 

Set  EC  (equation  count)  ■  0 

Initialize  putavay 

Skip  modifiers 

Increment  EC  by  1 

Increment  putavay 

Reset  ho 

Clear  TS 

Shift  EC 

Extract  into  b  pickup 

Jump  modifiers 

Increment  h 

Decrement  b 

Load  h 

Load  b 

Multiply 

Load  TS 

Add 

Putavay  summation 

If  b  greater  than  0,  repeat 

No,  summation  to  ansver  storage 

If  (n  ♦  r)  greater  EC,  0256 

Clear  tally,  Inversion  completed 


1   0 
1 


I'j«  in  M3  storage 

CC  storage 

I*   in  M1  storage 

M  storage 

Testvord  for  0131 

RS  storage 

Shift  Control 

Clear  Routine  (2001) 

M3  extractor 

Clear  Routine  (2003) 

Constant  for  01l6 

Constant  for  Clear  Routine 

CM  storage 

Constant  for  0113  etc. 


Ill  -  1  -  3 


MATRIX  INVERSION  PROGRAM   (CRC  -  102A) 
ADD   INST    M1      M2      M3  REMARKS 


0340 

00 

0000 

0000 

0000 

Available 

031*1 

00 

1000 

0000 

0000 

Constant  for  0156 

0342 

10 

0100 

0000 

0000 

Not  required 

03^3 

00 

0000 

0000 

0004 

Not  required 

03  44 

00 

0000 

0000 

0004 

Not  required 

0345 

00 

0000 

0000 

0004 

Not  required 

03  46 

30 

0364 

2100 

0000 

Reset  for  0152 

0547 

36 

2100 

2000 

0000 

Reset  for  0244 

0350 

00 

0000 

2100 

2001 

Testword  for  0275 

0351 

00 

1700 

0000 

1000 

Putavay  reset 

0352 

00 

0000 

0000 

0001 

EC  tally 

0353 

34 

3000 

2100 

0161 

Modifier  for  0252 

035^ 

34 

3000 

2100 

0251 

Modifier  for  0247 

0355 

00 

0003 

0000 

0000 

N  in  M1 

0356 

00 

0511 

0000 

0000 

J^ally 

0357 

00 

0000 

0000 

0002 

0360 

00 

0001 

0000 

0000 

Constant  for  0146 

0361 

00 

0000 

0000 

0030 

Shift  control  word 

036s 

00 

0500 

0000 

0000 

Constant  for  0126 

0363 

00 

7777 

0000 

0000 

M1  extractor 

03  64 

02 

2314 

1540 

0207 

Temporary  storage  (TS) 

0365 

00 

0000 

0000 

1702 

S  storage 

0366 

00 

0000 

0000 

0001 

(n  •♦•  r)  storage 

0367 

00 

4000 

0000 

0001 

"1"  in  PackBinFltPt 

0370 

26 

2100 

2100 

0500 

Clear  Routine  (2000) 

0371 

00 

0000 

0000 

0003 

N  storage 

0372 

34 

2004 

2000 

2000 

Clear  Routine  (2002) 

0373 

00 

0000 

0000 

0003 

RC  tally 

0374 

00 

2100 

2100 

1400 

Clear  Routine  (2004) 

0375 

00 

0077 

2100 

1777 

Testword  for  0105 

0376 

00 

0001 

0000 

0001 

Constant  for  0104 

0377 

30 

0000 

2100 

1700 

Not  required 

Channel  4  contains 

packed  floating  point 

arithmetic  routines 

0400 

30 

3000 

0443 

2006 

Prepare 

0401 

32 

2006 

0475 

0433 

exit 

0402 

30 

2001 

2100 

2007 

Store  2nd  operand 

0403 

52 

2000 

0434 

2102 

Extract  sign  of  exp1 

0404 

27 

2002 

0451 

2002 

Shift  to  operating  location 

0405 

32 

2000 

0453 

2002 

Extract  mag  of  exp1 

0406 

52 

2000 

0444 

2103 

Extract  sign  and  mag  of  mag1 

0407 

31 

2002 

2003 

2000 

Shift  to  2000,2001 

o4io 

32 

2007 

0434 

2102 

Extract  sign  of  exp2 

04ll 

27 

2002 

0451 

2002 

Shift  to  operating  location 

04l2 

52 

2007 

0436 

2002 

Extract  mag  of  exp2 

04l5 

32 

2007 

0444 

2103 

Extract  sign  and  mag  of  mag2 

04l4 

36 

2006 

0435 

2006 

Pickup  entry  command 
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MATRIX  INVERSION 

PROGRAM 

(CRC  -  102A) 

ADD 

INST 

M1 

M2 

M3 

REMARKS 

04l5 

30 

2006 

0437 

2006 

Shift  to  M1 

04l6 

32 

2006 

046l 

0417 

Set  0417 

0417 

32 

0273 

046l 

2107 

Extract  control 

0420 

54 

2007 

0400 

0440 

If  control  -^  5000,  0440 

0U21 

33 

2002 

2000 

0467 

No,  operations 

0U22 

36 

2002 

2000 

2006 

for  normal 

0423 

30 

2003 

2006 

2007 

floating  point 

0424 

55 

2001 

2007 

2001 

addition 

0425 

37 

2001 

3000 

0462 

routine 

0426 

31 

2000 

2001 

2002 

Store  in  2002,  2003 

0427 

27 

2002 

0452 

2104 

Shift  sign  of  exp  for  pack 

0U30 

32 

2002 

0436 

2004 

Extract  mag  of  exp  for  pack 

0U3I 

32 

2004 

0476 

2003 

Pack  exponent 

OU32 

30 

2003 

2100 

2000 

Packed  result  to  2000 

OU35 

34 

3000 

2100 

0274 

Exit  to  main  program 

0U3U 

00 

0000 

0000 

0100 

Extractor  for  sign  of  exp 

0435 

00 

0000 

0000 

0001 

Constant  for  04l4 

OU36 

00 

0000 

0000 

0077 

Extractor  for  mag  of  exp 

OU37 

00 

0000 

0000 

0030 

Shift  control  word 

o44o 

34 

2007 

0450 

0445 

If  control  ^  3200  0445 
Change  sign  of  2   operand 

044l 

56 

2100 

2003 

2003 

0442 

34 

3000 

2100 

0421 

Transfer  to  addition 

0U3 

02 

0000 

0000 

0030 

Shift  control  word 

0444 

02 

7777 

7777 

7600 

Extractor  for  magnitude 

0445 

34 

2007 

0474 

0454 

If  control  ^  3400,  0454 

0446 

35 

2000 

2002 

2000 

No,  operation  for  normal 

0447 

25 

2001 

2003 

2001 

floating  point 

0U50 

34 

3200 

2100 

0426 

multiplication 

OU51 

00 

0000 

0000 

0037 

Shift  to  pack  commands 

0452 

02 

0000 

0000 

0037 

Shifter  for  sign  of  exp 

0453 

00 

0000 

0000 

0077 

Extractor  for  mag  of  exp 

0U5U 

36 

2000 

2002 

2000 

Operations  for 

0455 

23 

2001 

2003 

2001 

normal  floating  point 

0U56 

37 

2001 

3000 

0462 

division 

OU57 

54 

3000 

2100 

0426 

Transfer  to  pack  commands 

0460 

02 

0000 

0000 

0001 

Constant  for  0462 

046l 

00 

3777 

0000 

0000 

M1  extracotor 

0462 

27 

2001 

0460 

2001 

Normal 

0^63 

30 

2001 

0460 

2001 

correct 

0464 

27 

2001 

0477 

2001 

overflow 

0465 

35 

2000 

0477 

2000 

procedure 

0466 

34 

5000 

2100 

0426 

Transfer  to  pack  commands 

0467 

56 

2000 

2002 

2006 

Operations  for 

0470 

30 

2001 

2006 

2007 

normal  floating 

0471 

30 

2002 

2100 

2000 

point  addition 

0472 

35 

2003 

2007 

2001 

when  2002  >  2000 

0473 

37 

2001 

3000 

0462 

see  0421 

0474 

34 

3400 

2100 

0426 

Transfer  to  pack  commands 

0475 

00 

0000 

0000 

7777 

M3  extractor 

0476 

00 

0000 

0000 

0177 

Sign  and  exponent  extractor 

0477 

00 

0000 

0000 

0001 

Constant 
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FLOW  CHART  OF  MATRIX   INVERSION  AND  SOLUTION 


Add  N+l   to 
"1"  PA 

Cotanand 


T 


<*.  j =*Add   1   to  CM 


©- 


(,      "Matrix"       [H]  J 


X 


Clear  0500-1777 

(H|— >■  1700 

Set  RC  =  0 

Set  N+l  in  modifier 

Initialize  "1"  PA 


KD 


This  N  is  the 
order  of  the 
matrix,  viz.  n 


'l"-»a 


M, 


'(i-J) 
add  1  to  RC 


Yes 


Unit  Matrix 


Set  1->CM 
1701  +  n  +  r->S 


■CZiZZZ) 


Reset  PA  to  0000 
Set  RC  to  0 
Reset  a   to  a.,  pickup 


Setup  to 
multiply 
in  rows 


S  -  CM-*RS 
Clear  CC 
Clear  TS 


Row  selection 


Yes 


Yes 


No 


[Aj  •   [RS>  B 

prs]  -»■  Qao-»TS 


Add   1   to   CC 


No 


Putaway  ^TS  ->  0000 
(initially) 
Add   1  to  RC 


Yes 


Yes 


Add   1   to  a 


iJ 


pickup 
Add   1   to  putaway 


[cm]  -l-»  J 
(jJ'N  +  0500  -  IL 
i 
0000  +  [jj  =  P. 


Set   CC   to  0 


r  r  r 

Add    1   to  CC 


CM  >  Cc\_X?i_ 


-> 


RS   +    l->  RS 
Add    I    to  a        pickup 


Mod   to   I£      +  Pk_ >  I£ 
r  r  r 


> 


r 


r 


Set  M  to  0 


Yes 


0500  +  M-*I. 

J 

Set  P. — »  0000 
k 

r 

Set  Dec  to  N+l 
Set  RC  to  1 


Yes 


0 


No 


^ 


Add  1  to  RC 


Yes 


v\ 


Inc  I.  to  I- 
r     r 

Add  1  to  M 


No 


add  1  to  EC 
Mod  PA 


Add  1  to  CM-4»CM 
Reset  a   pickup  to  a.. 

PA  "1"  NPFP-*-  0000 

Set  PA  to  0001 
Alter  "out"  cmd  to  rtn  to  8 


Restore  "out"  cmd 
Set  EC  to  0 
Set  putaway 


Reset  hE(, 

Clear  TS 

Shift  EC 

Initialize  b  pickup 


TS  +  (hEC-i),(bi)-^TS 


i  +  1 


No 


put  away 


Yes 


TW  .  000021002001 


Yes 


END 


III  -  2  -  1 


APPENDIX  IV 
LAG RANG IAN  INTERPOLATION  METHOD  FOR 
THE  SOLUTION  OF  ALGEBRAIC  EQUATIONS1 

The  Lagrangian  interpolation  method  described  by  D.  E.  Muller  was 
selected  for  the  solution  of  large  degree  polynominals .   This  procedure 
was  programmed,  with  certain  modifications,  and  tested  on  the  CRC-10ZA 
computer.   By  removing  the  roots,  as  found,  by  synthetic  division,  rather 
than  by  the  algorithm  recommended  by  Muller,  accuracies  of  at  least  eight 
decimal  places  were  obtained  for  all  roots  of  equations  up  to  degree  21. 
This  was  considerably  better  than  the  accuracies  reported  by  Muller,  and 

incidentally,  probably  entailed  considerably  more  iteration. 

-13 
A  graph  of  solution  time  vs.  degree  of  equation,  for  (■ ■_  *     10    ,  is 

presented  in  Figure  IV-1.   Convergence  to  a  root  generally  requires  be- 

-13 
tween  7  and  14  iterations  (for  the  £  of  10    ).   The  iterative  loop  is 

completed  in  about  1.5  minutes,  varying  slightly  (not  more  than  10%) 

according  to  the  degree  of  the  equation.   In  making  these  tests,  it  was 

found  that  when  equations  had  roots  of  the  same  modulus,  there  is  a  pro- 

24 
nounced  slow  down  in  the  rate  of  convergence.   The  equations  x   +  1  ■  0, 

20 
and  x   4-1*0  required  about  12.5  minutes  per  root,  whereas  the  more 

general  polynominal,  (equation  A  under  Results  of  Tests)  only  required 

about  9.4  minutes  per  root.   Round-off  error  was  not  present  in  the 

x   +  1  equation  and  affected  only  two  roots  in  the  solution  of  x   +1*0 

within  the  eighth  decimal  place. 

D.  E.  Muller,  Solving  Algebraic  Equation  by  an  Automatic  Computer, 
Mathematical  Tables  and  Other  Aids  to  Computation,  Vol.  X,  No.  56,  Oct., 
1956,  pp.  208-215. 
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Thus  far  the  program  has  failed  only  once.   When  solving  an  equation 
with  a  real  double  negative  root  and  a  pair  of  complex  roots  of  about 
the  same  modulus,  viz.,  4;  the  rate  of  convergence  was  normal  until  about 
four  decimal  places  had  been  established  and  decreased  so  as  to  become 
barely  perceptible.   The  program  did,  however,  remove  all  roots  interior 
to  modulo  4  without  difficulty. 

The  basic  outline  of  the  Lagrangian  interpolation  method  is  as 
follows : 

Successive  iterations  toward  any  particular  root  are  obtained  by 
finding  the  nearer  root  of  a  quadratic  whose  curve  passes  through  the 
points;  z^  f^);  z^,  fCz^);  zi.2»  f(zi-2^' 

The  Lagrangian  formula  may  be  written  as  L.  I  f(z)J  ■  b_z  +  b.z  +  b  ; 
where  the  b's  satisfy  the  system  of  equations: 

Vi2  +  bl2i  +  b2  '  f(zi} 

Vi-1  +  blZi-l  +  b2  "  f(zi-l} 

Vi-2  +  Vi-2  +  b2  "  f(zi-2> 
and  these  equations,  of  course,  become  almost  one  in  the  same  as  the 
iterative  process  converges,  depending  on  the  accuracy  which  we  specify. 

We  are  to  construct  a  polynominal  (quadratic)  whose  curve  passes 

2 
through  the  three  points.   By  the  Lagrangian  interpolation  formula  ,  we 

write: 

r       "l                 (*  •  «j   i)(z  "  z<_o)                           (z  "  z<)(z^"z<   ?) 
L,      f(z)   -  f(zj  t w « r-  +   f(z4    .) * \.    w - 

iL       J  ^    (zi    ■   Zi-l)(zi    '   Zi-2>  i"1      (zi-l    *   Zi)(zi-1    "   zi-2> 

(z    -   zi)(z    -   zi.1) 

+   f(V2>   <\-2   ■  zi><zt-2    -  *i-l> 

2 
W.  E.  Milne,  Numerical  Calculus,  Princeton  University  Press,  Princeton, 

N.  J.,  1949,  for  example. 
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If  we  make  the  following  simplifying  substitutions: 


z    -  z     ■  h 


Zi   "  Zi-1  '  hi 


Zi-1   "   Zi-2   *  hi-l 


h/h.    «    X 

i 

Vhi-i  m\ 

1  +\i  -  Ai 


we  may  express  the  coefficient  of  f(z.)  as: 

since        (z    -  z1mI)(«   -   zim2)     =     ^  +   ^    (z    -   z._2) 
(Z1   *  Zi-l)(zi   "Zi-2)     "  (zi   '  Zi-2> 


1  Zi-1   "  Zi-2 


Ai  i 


,   and 


i-1 


+  X  +  i 


z    -   z 


\ 


4i 
Xi 


^ 


X 


*=»-  -^  ♦  X 


zl  -  Vl     Xi 


z    -   z 


i-2 


Zi-1    "   Zi-2 


s%[fr\ 


z    -   z 


i-2 


Zi   '   Zi-2 


and   finally 


;    therefore, 


\  \\  i  (z    "   zi-i)(z    "  zi-?> 

(1  +  X)  (1  AX. a;1) 


t- — — rr- — — r     the   coefficient  of   f(z    ), 

Ui        Zi-1A    i  i-2;  1 


which  when  expanded   is: 


$  Xi  a:1  +xXi  a^1  +  x- 1  -  x2  \t  a^1  ♦  Xa:1^  -ap  + 1 


(QED) 


By  similar  manipulation,  the  coefficients  of  f(z    )  and  f(z   .)  cart  be 
found  and  the  desired  quadratic  in  A  is: 
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\  [«(o]     -  X2  A*1 


f(zi-2}    Ai    "  f(z 


t-1)     \\      ♦   fCz^^ 


Xa_1  [f(^i_2)Xi2  -  f<zi-i>  Ai2  +  fzi^i  +Ai)J  +  f(zt> 

A  single  iterative  step  is  obtained  by  letting  z    be  such  a  value 
of  2  that  L,  U(z)   *  0«   We  may  solve  the  quadratic  inversely  to  obtain 

X  which  is  A 


Zi+1  "  2i 
1+1  "  Zi  *  Zi-1 


from  which  we  obtain  the  new  z   and  repeat 


the  iterative  process. 

All  operations  are  in  complex  numbers  and  in  floating  point  arithmetic, 
A  root  is  taken  when 


z,  ,  -  z . 
i+1    i 


<  €. 


For  proof  of  convergence  of  the  process  the  reader  is  referred  to  Dr. 
Muller's  paper. 
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RESULTS   OF  TESTS 


Equation  A 
16 


-  9.0  z15  +  34.75   z14 


13  12  11 

102.25   z       +  264.0  z        -  479.0  z 


+   771.5   z10   -   1154.5   z9  +  978.0  z8   -   1225.0  z?  +  912.75   z6 

+  254.75   z5  +   1411.0  z4  +   1319.0  z3  +  1022.5   z2  +  470.5  z 

-44  -13 

+   105.0      6-2        (=   10        ) 

Total  solution  time:   2  hours,  34  minutes,  10  seconds  on  CRC-102A. 


Roots 


-0.284460054 

+ 

J 

0.598759363 

-0.376766336 

+ 

J 

0.188840569 

-0.1012  79181 

+ 

J 

1.429775042 

-0.044975589 

+ 

J 

1.738162299 

+1.999999996 

-0.011258189 

+ 

J 

1.933289317 

+2.500000009 

-0.181260648 

+ 

J 

1.034830851 

+3.000000001 

+3.499999991 

IV 


RESULTS  OF  TEST  (Continued) 


Equation  B 

30 

x   + 

1-0 

€   -  2-44  (-  io'13) 

Total 

solution  time: 

6  hours,  23  minutes,  30  s 

Roots : 

-0.000000000 

+ 

)     0.999999999 

-0.951056516 

+ 

)     0.309016994 

+0.866025403 

+ 

j     0.500000000 

-0.587785252 

+ 

j     0.809016994 

+0.994521895 

+ 

j     0.104528463 

-0.207911690 

+ 

j     0.978147601 

-0.994521895 

+ 

j     0.104528463 

+0.406736637 

+ 

j     0.913545463 

+0.207911688 

+ 

|     0.978147598 

+0.743144829 

+ 

1     0.669130603 

-0.743144826 

+  : 

|     0.669130607 

+0.587785258 

+ 

j     0.809017000 

-0.866025404 

+  j 

|     0.499999998 

+0.951056517 

+  : 

1     0.309016993 

-0.406736643 

+  * 

|     0.913545458 
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IV 


LAGRANGIAN  INTERPOLATION  METHOD  FOR  HIGH  DEGREE  POLYNOMIALS 


ADD   INSI 


MJ 


M2 


M3 


REMARKS 


Channels  0  and  1  are  used  temporary  storage. 

Channel  2  contains  a  decimal  to  binary  floating  point  conversion 
which  while  not  needed  in  overall  program  is  included  here. 


0200 

26 

15^0 

0250 

2005 

Build  and 

0201 

55 

0246 

2005 

0255 

store  testword 

0202 

55 

0251 

2100 

0210 

Reset  pickup 

0205 

55 

0255 

2100 

0211 

Reset  pickup 

020U 

52 

2100 

0252 

0270 

Reset  putaway 

0205 

05 

5000 

5000 

5000 

Clear  buffer 

0206 

50 

025^ 

2100 

0265 

Reset  extractor 

0207 

52 

0251* 

0255 

2004 

Set  2001*  to  1 

0210 

50 

0016 

2100 

0250 

Pickup  integer 

0211 

50 

0017 

2100 

021*7 

Pickup  fraction 

0212 

5^ 

0250 

2100 

02  lU 

Is  Integer  *  0? 

0215 

5U 

5000 

2100 

0251 

No,  test  fraction 

021^ 

52 

0250 

0255 

0227 

Yes,  extract  sign 

0215 

52 

0250 

0265 

2005 

Extract  least  slg.  digit 

0216 

26 

2001* 

2005 

2006 

L.S.D.  x  (2001*) 

0217 

55 

2006 

2001 

2001 

Accumulate  products 

0220 

26 

2001* 

0257 

2001* 

Increase  2001*  by  power  of  10 

0221 

50 

0250 

0262 

0250 

Shift  decimal  word  1*  right 

0222 

3h 

0250 

2100 

0215 

If  conversion  not  complete,  0215 

0225 

51 

2000 

2001 

2000 

If  complete  normalize 

0221* 

55 

2000 

0261 

2000 

and  correct  exponent 

0225 

5U 

021*7 

2100 

0255 

Is  fraction  ^  0? 

0226 

5* 

5000 

2100 

0267 

No,  skip  to  putaway 

0227 

00 

0000 

0000 

0000 

Sign  storage 

0250 

00 

0000 

0000 

0000 

Integer  storage 

0251 

5^ 

02 1*7 

2100 

0251* 

I  *  0,  Is  fraction  ^  0? 

0252 

5* 

5000 

2100 

0267 

If  both  are  0,  putaway  0 

0255 

00 

0016 

2100 

0250 

Testword 

025** 

52 

021*7 

0255 

0227 

Extract  sign  when  fraction  only 

0255 

52 

021*7 

0265 

2005 

Extract  least  slg.  frac.  dec.  dig 

0256 

25 

2005 

0265 

2005 

Divide  by  10 

0257 

50 

0265 

0256 

0265 

Shift  extractor 

02l*0 

52 

021*7 

0265 

2005 

Extract  next  digit 

021*1 

25 

2005 

0265 

2005 

Divide  by  10 

021*2 

54 

0260 

0265 

0257 

If  conversion  not  complete,  0257 

02h  5 

5^ 

2001 

2100 

0266 

If  integer  x  0,  transfer  to  add 

021*1* 

51 

2002 

2005 

2000 

If  integer  was  0,  normalize 

021*5 

51* 

5000 

2100 

0267 

and  putaway 

021*6 

00 

0002 

2100 

0250 

Testword  addend 

02U7 

00 

0000 

0000 

0000 

Fraction  storage 

0250 

00 

0002 

0000 

0000 

Testword  build  constant 

0251 

50 

0000 

2100 

0250 

Integer  pickup  reset 

0252 

00 

0000 

0000 

7777 

M3  extractor 
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LAGRANGIAN  INTERPOLATION  METHOD  FOR  RICH  DEGREE  POLYNOMIALS 
ADD   INST    M1      M2      M3  REMARKS 


0255 

50 

0001 

2100 

02U7 

Fraction  pickup  reset 

025H 

00 

0000 

0000 

0017 

Extractor  reset 

0255 

02 

0000 

0000 

0001 

Sign  extractor  and  -1 

0256 

00 

0000 

0000 

ooou 

Shift  control  -  U  left 

0257 

00 

0000 

0000 

0012 

"Power  of  10" 

0260 

00 

7fcOO 

0000 

0000 

Extractor  testword 

0261 

00 

0000 

0000 

oouu 

To  correct  exponent 

0262 

02 

0000 

0000 

ooou 

Shift  control  -  k   right 

0263 

00 

5000 

0000 

0000 

10  for  division 

026*4 

00 

0000 

0000 

0002 

2 

0265 

00 

7U00 

0000 

0000 

Extractor  working  storage 

0266 

5^ 

3000 

2100 

15  Uo 

Transfer  to  addition 

0267 

52 

0227 

0255 

2001 

Extract  sign  into  mag  portion 

0270 

50 

2000 

2100 

0016 

Putaway  exponent 

0271 

56 

0270 

0255 

2007 

Build  magnitude 

0272 

52 

2007 

0252 

0275 

putaway 

0275 

50 

2001 

2100 

0015 

Putaway  magnitude 

0271* 

55 

0210 

0250 

0210 

Modify  integer  pickup 

0275 

55 

0211 

0250 

0211 

Modify  fraction  pickup 

0276 

55 

0270 

0261* 

0270 

Modify  putaway 

0277 

5^ 

0255 

0210 

0205 

Have  all  numbers  been  converted? 

Channels  3  through  7,  and  1000  through  1012  contain  main  program 
for  Lagrangian  Interpolation  Method;  in  the  program  detailed 
here  Channels  10  and  11  contain  a  binary  floating  point  to 
decimal  conversion  routine  which  would  not  be  required  for 
the  overall  program.  Channel  13  is  used  for  tempoary  storage 
and  Channels  Ik   through  16  contain  the  floating  point  arith- 
metic routines  for  complex  numbers „ 
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LAGRANGIAN  INTERPOLATION  METHOD  FOP  HIGH  DEGREE  POLYNOMIALS 


A  DD 

INST 

Mx 

M2 

M3 

0300 

50 

13^0 

2100 

1372 

0301 

05 

3000 

5000 

0000 

0302 

04 

3000 

3000 

0100 

0503 

35 

0301 

1376 

0301 

0304 

35 

0302 

1376 

0302 

0505 

34 

137^ 

0302 

0301 

0306 

56 

0301 

1^57 

0301 

0307 

36 

0302 

1437 

0302 

0310 

30 

1531 

2100 

2000 

0311 

36 

2100 

1337 

2001 

0312 

31 

2000 

2001 

1300 

0313 

31 

1333 

1337 

1320 

03lfc 

31 

1531 

1357 

1301 

0*>i5 

31 

1333 

1337 

1321 

0316 

31 

1333 

1337 

1302 

0317 

31 

1335 

1357 

L322 

0320 

55 

2100 

2100 

1303 

0521 

36 

2100 

1337 

13U 

0322 

31 

1333 

1337 

1325 

0323 

26 

1372 

l466 

2007 

0324 

26 

2100 

2100 

1434 

0325 

26 

2100 

2100 

1455 

0526 

35 

06 1*0 

2007 

0327 

0327 

31 

0102 

0103 

2000 

0330 

31 

2000 

2001 

1307 

0331 

31 

1333 

1337 

1327 

0332 

56 

0327 

1473 

200? 

0333 

32 

2007 

1464 

0554 

0334 

31 

0070 

0071 

2002 

0335 

34 

3000 

2100 

15^0 

0356 

36 

033** 

1473 

0334 

0357 

5^ 

0640 

033^ 

03^1 

0340 

34 

3000 

2100 

0334 

034l 

34 

1475 

2100 

O3U7 

0342 

31 

2000 

2001 

1^34 

03^3 

36 

0327 

1U66 

03  44 

0544 

31 

0100 

0101 

2000 

0345 

36 

03  44 

1473 

2007 

0346 

3^ 

3000 

2100 

0333 

0347 

31 

2000 

2001 

1^35 

0350 

31 

lU5i* 

1475 

2002 

0351 

34 

3000 

2100 

1540 

0352 

31 

2000 

2001 

1306 

0353 

31 

1333 

1357 

1326 

0354 

31 

1454 

1^75 

2000 

0355 

31 

1435 

1U76 

2002 

0356 

34 

3000 

2100 

1516 

0357 

31 

2000 

2001 

1305 

REMARKS 


N  to  N' 
Transfer 
the  contents 
of  Channel  0 

to 
Channel  1 

and 
reset. 

Initial  setup 
for  z^,2 

Initial  setup 

for  Zj„i 
Initial  setup 

for  z± 
Initial  setup 

for  Xf 

-  1/2  +  JO 
Commence  pickup  build 
Clear  1434,  1475 
Clear  1455,  1476 
Set  pickup  for  A 
Pickup  A0 
Ao  is  fkof  f(zi) 
&*f(zi)  «  0 
Build  A0  pickup 

and  enter  in  0354 
An  pickup  command 
Transfer  to  floating  addition 
Modify  An  pickup  for  A^g 
If  sum  complete,  putavay 
No,  repeat  to  0334 
Has  even  sum  putavay  been  made? 
No,  even  sum  to  1434,  1475 
Build  A^^  pickup 
Ar-1  pickup 
Set  An  pickup  for  odd  terms 

and  repeat  back  for  odd  sum 
Putavay  odd  sum  in  1435,  l476 
Load  even  sum 
Odd  *  even  for 

f(«!.l) 

and  clear  imaginary 
Load  even  sum 
and  odd  sum 
for  subtraction 
to  form  f (zj.2^ 


IV  -  1 


LAGRANGIAN  INTERPOLATION  METHOD  FOR  HIGH  DEGREE  POLYNOMIALS 


ADD   INST 


MJ 


M2 


M3 


REMARKS 


0360 
O361 
0362 

0563 
O36U 
O365 
O366 
0367 
0370 
0371 
0372 
0573 
037^ 
0375 
0376 
0377 

01*00 

Ol*01 
0l*02 
0^03 
Ol+OU 
0^05 
0k06 
0^07 
0*10 

oim 
01*12 

Ol*13 

0U1U 

OlH5 

01*16 
01*17 
0U20 

Ol*21 
01*22 
01*23 
0U21* 
01*25 
01*26 
01*27 
01*30 
01*31 
01*32 
01*33 
Ol*3l* 
01*35 
01*56 
OU37 


31 
51 
31 
31 
31 
3* 
31 
31 
31 
31 
3^ 
01* 

31 
31 
31 
31 
31* 
01* 

31 
31 
^ 
31 
31 
3^ 
51 
31 
3^ 
31 
31 
3I* 

31 
31 
35 
35 
3* 
Ol* 

05 
31 
31 
3* 
01* 

05 
31 
31 
3U 
Ol* 
31 
31 


1333 
1303 
1323 
1331 
1335 
3000 
2000 
2002 
1306 
1326 
3000 
3000 
1305 
1325 
1303 
1523 
5000 
3000 
1307 
1527 
3000 
1330 

1332 

3000 

130U 

1321* 

3000 

1303 

1323 

5000 

1307 

1327 

2001* 

2006 

3000 

3000 

3000 

1301* 

1321* 

3000 

3000 

3000 

1303 

1323 

3000 

5000 

1303 

1323 


1337 
13W 
1361* 

1335 
1337 
2100 
2001 
2003 
131*7 
1367 
2100 
3000 
131*6 
1366 
13l*l* 
1361* 
2100 
3000 
1350 
1570 
2100 
1371 
1375 
2100 

13^5 

1565 

2100 

13l*l* 

1361* 

2100 

1350 

1370 

0261* 

0261* 

2100 

3000 

3000 

13U5 

1365 

2100 

3000 

3000 

131*1* 

1361* 

2100 

3000 

13l*l* 

1361* 


1325 
2000 
2002 
2001* 
2006 
ll*10 
1301* 
132U 
2001* 
2006 
1607 
1330 
2000 
2002 
200U 
2006 
1607 
1310 
2001* 
2006 

1U10 
2001* 
2006 

11*16 

2001* 
2006 
1607 
200U 
2006 
1607 

2001* 
2006 
2001* 
2006 
1607 
ll*l*0 
1330 
2001* 
2006 
1607 
1330 
1310 
2001* 
2006 
1607 
1310 
2000 
2002 


and  clear  imaginary 
Load 

Xi 
Load 

1  +  JO 

Transfer  to  cpx.   add. 

Put away 

Load 

f(zi-l)' 

Transfer  to  cpx.  mpy. 

Store  in  B 

Load 

f(zi-2) 
Load 

>1 

Transfer  to  cpx.  mpy. 

Store  in  A 
Load 

f (z,   2)  \i  +  f(z1) 
Load 

Lf(zi-2Mi  -  f(zi^1)^i  +  ffzjU 
Load 

Load 

ffzl-Aj-  At 
Load 

t(%±) 

and   "multiply 

by  1* 
Uf(z1)A1Xi  f(Z)«  U*AC 
Store  in  C 
Load  f  fzlral)A  i 
Load 

^i  <L 

f(«l.l)-&l 

Store   In  B 

f fzi.g)  X  i  from  A 

Load    Xi 

Multiply  for  f(Zi„2)  Xj2 

Store  in  A 

Load 


TV 


LAGRANGIAN  INTERPOLATION  METHOD  FOR  HIGH  DEGREE  POLYNOMIALS 
ADD   INST    M1      M2      M3  REMARKS 


0440 
044l 
0442 
0443 
0444 
0445 
0446 
Okkj 
0450 
OVjl 
0^52 

0453 
0454 

0455 
OU56 

OU57 
0460 
046l 
0462 
0U65 
0464 
0465 
0466 
0467 
0U70 
0471 
OU72 

0473 
0474 

0475 
0476 

0477 
0500 
0501 
0502 
0505 
0504 
0505 
0506 
0507 
0510 
0511 
0512 
0513 
0514 

0515 
0516 

0517 

0520 


31 
31 
34 
31 
31 
3* 
31 
31 
54 
31 
31 
34 
04 

31 
31 
34 
31 
31 
34 
04 

31 
30 
36 
3* 
34 
31 
34 
31 
31 
31 
34 
31 
31 
3* 
56 
34 
31 
34 
31 
31 
31 
34 
36 
34 
31 
34 
31 
30 
30 


1304 

1323 
3000 
1307 
1327 
3000 
1310 
1312 
3000 
1330 
1332 
3000 
3000 
2000 
2002 
3000 
1440 
1442 
3000 
3000 
2000 
2002 
2100 
3000 
3000 
2000 
3000 
1434 
2000 
1440 
3000 
2000 

1331 
3000 
2000 
3000 
1454 
3000 
2000 
1455 
1331 
3000 
2000 
3000 
1434 
3000 
2000 
2000 
2000 


1345 
1365 
2100 
1350 
1370 
2100 
1351 
1353 
2100 
1371 
1373 
2100 
3000 
2001 
2003 
2100 
1401 
1403 
2100 
3000 
2001 
2100 
2003 
2100 
2100 
2001 
2100 

1475 
2001 
1401 
2100 
2001 
1337 
2100 

1331 
2100 

1475 
2100 
2001 
1476 

1337 
2100 

1331 
2100 

1475 
2100 
2001 
2100 
2100 


2004 
2006 
1410 
2004 
2006 
1607 
2004 
2006 
1410 
2004 
2006 
1416 
1330 
2004 
2006 
1607 
2004 
2006 
1416 
1440 
2004 
2006 
2007 
1607 
1551 
1434 

1551 
2004 
1434 
2000 
1507 
3*35 
2002 
1540 
2000 

1551 
2004 
1500 
1436 
2002 
2000 
1516 
2000 
1551 
2004 
1500 
1434 
1451 
1453 


Load 

Load 
ffzi) 

f(zi)    (  A1   +  At) 
Load 

2 


+v 


f(Zi-2>XiS 


f(Zi-2^V~  f^zi-l)A12  +  f(ziKAi^i) 


Result  ■  g±;   store  in  B 
Load  B  for  complex 

multiplication,   i.e.   B2 
gi2     OR  B2 
Load  4AC 
B2  -  4AC 
Store  in  C 
Load 

B2  -  4AC 
for  conjugate 
multuplication 
(B2  -4ACHB2  -  4AC)*  =  r2 
Compute  r 
Store  r  in  TSX 
Compute  -J"?"" 
Load  r 
Store  -jr*    in  TS 

(B2  -  4AC)  =  x 
x/r  ■  cos  8 
Store  cos  6   in  TS2 
Load  1 
1  +  cos  9 
(1  +  cos  9  )/2 
-|(1  *  cos  e)72     -  cos  0/2 
Load  ^r 

Jr   (cos  9/2)  =   *• 
*'  to  TS3 
Load  cos  6 
Load  1 
1  -  cos  & 
(1  -  cos0)/2 

{(1  -  cos  9)/2     *  sin  0/2 
Load  -JF 

■Jr    (sin  0/2)  -   a1 
A*  to  TSi  r 

Putavay(2>0 

and    /3X  exps 


IV  -  1 


LAGRANOIAN  INTERPOLATION  METHOD  FOR  HIGH  DEGREE  POLYNOMIALS 
ADD   INST    M1      M2      M3  REMARKS 


0521 

0522 

0523 

052* 

0525 

0526 

0527 

0530 

0531 

0532 

0533 

05  5* 

0535 

0536 

0537 

05*0 

05*1 

05*2 

05*3 

05** 

05*5 

05*6 

05*7 

0550 

0551 

0552 

0553 

055* 

0555 

0556 

0557 

0560 

0561 

0562 

0563 

056U 

0565 
0566 
0567 
0570 
0571 
0572 
0573 
057* 

0575 
0576 

0577 


30 

30 

52 

50 

36 

32 

33 

50 

36 

3* 

30 

36 

05 

31 

31 

3* 

31 

30 

36 

3* 

31 

05 

31 

31 

3* 

31 

30 

36 

5* 

51 

05 

33 

33 

3* 

31 

31 

3* 

31 

31 

3* 

0* 

31 
31 
35 
36 

35 
36 


1*36 
1*36 
1*75 
2007 
2100 
1*77 
1*03 
2007 
2100 
5000 
2007 
2100 
3000 
1*50 

1*51 
5000 
2000 
2002 
2100 
5000 
2000 
3000 
1*52 
1*53 
3000 
2000 
2002 
2100 
3000 
2000 
3000 

1*3* 
1*36 
1*75 
1*51 
1*50 
3000 

1*53 
1*52 

3000 
3000 
1307 
1527 
130* 
2100 
132* 
2100 


2100 
2100 

1575 

2100 

2007 

1575 

2100 

2100 

2007 

2100 

2100 

2007 

3000 

1*60 

l*6l 

2100 

2001 

2100 

2003 

2100 

2001 

3000 

1*62 

1*63 

2100 

2001 

2100 

2003 

2100 

2001 

5000 

1*36 

1*3* 

1*77 

l*6l 

1*60 

2100 

1*63 

1*62 

2100 

3000 

1350 

1370 

1331 

13*5 

1331 

1365 


1*50 
1*52 
2107 
1*61 
1*63 
2107 
0533 
1*62 
1*60 

0535 

1*60 

1*62 

1330 

200* 

2006 

1*16 

200* 

2006 

2007 

1607 

1*3* 

1330 

200* 

2006 

1*16 

200* 

2006 

2007 

1607 

1*36 

1330 

0566 

O563 

0566 

2006 

200* 

0570 

2006 

200* 

1*10 

1330 

2000 

2002 

200* 

2005 

2006 

2007 


Put  away  o{Q 

and  =\x  exps 
Set  /3C  mag  positive 

and 

/31  mag  negative 
Clear       mag  sign  bit 
If  JUB2  -  *AC)  *  0,  0533 
No,  putavay     <*i  positive 

and    °<0  negative 
Skip  to  0535 
Putavay  o<0  positive 

and     *!  negative 
Load  g1 
Load  <x0 
Load     (3© 
Subtract 
Load  for 

conjugate 

multiplication 
Multiply  and 

store  in  TSX 
Load  g± 
Load  <*! 
Load  (hi 
Subtract 
Load  for 

conjugate 

multiplication 
Multiply  and 

store   in  TS3 
Load  g 

If  TSj.  i»  TS3,  0566 
If  TS3  *  TSX,  O563 
If  TSj.  ^TS3,  0566 
Load  *<0 

and    fie 
or 


load 

B  +-J  B*   -*AC 
Store  in  B 
Load 

f(«±) 

Load 

-2  A 

for  1 

^A^fCzj) 


IV 


LAGRANGIAN   INTERPOLATION  METHOD  FOR  HIGH  DEGREE  POLYNOMIALS 


ADD       INST 


MJ 


M2 


M^ 


REMARKS 


0600 
0601 
0602 


3^ 
31 
31 


0603   3^ 
0601*  Ok 


O605 
O606 
0607 
0610 
06ll 
0612 
0613 
06ll* 
O615 
06l6 
0617 
0620 
0621 
0622 
0623 
0621* 
0625 
0626 

0627 
O63O 
0631 
0632 
0633 
O63I* 
0635 
0656 
0637 
061*0 
061*1 
061*2 

061*3 
061*1* 
061*5 
061*6 
061*7 
0650 
0651 
0652 
0653 
065I* 

0655 
0656 

0657 
0660 


31 
31 
31 
31 
3^ 
31 
31 
3*v 
31 
31 
3* 
Oh 

31 
31 
3^ 
31 
31 
^ 

33 
3U 
33 
26 

35 
30 
35 
31 
31 
31 
31 
32 
3* 
31 
32 

31 
3U 

31 
31 
31 
35 
3U 
31* 
31 
00 

01* 


3000 
1550 
1532 
5000 
3000 
1302 
1322 
1301 
1321 
3000 
1330 
1332 
3000 
1502 
1322 
3000 
3000 
1302 
1522 
3000 
1302 
1522 
3000 
Ik^h 
3000 
11*51* 
1372 
11*70 

ll*70 
2007 
1310 
1312 
0100 
1333 
0657 
3000 
2002 
0677 
0101* 
3000 
1602 
1310 
1312 
061*6 
061*6 
3000 
0101* 
0000 

3000 


2100 
1371 
1373 
2100 
3000 
13^3 
1363 
13l*2 
1362 
2100 
1371 
1573 
2100 

13^3 
1363 
2100 
5000 
13^3 
1563 
2100 
13^3 
1363 
2100 
2000 
2100 
2002 
11*66 
2007 
2100 
0261* 
1351 
1353 
0101 
1337 
1636 
2100 
2003 
1636 
0105 
2100 
161*3 
1551 
1353 
ll*66 
0656 
2100 
0105 
0000 
3000 


1607 
2001* 
2006 
1650 
1330 
2000 
2002 
2001* 
2006 
11*16 
2001* 
2006 
1607 
2001* 
2006 
ll*10 
1310 
2001* 
2006 
11*16 
2001* 
2006 
1650 
0631 
0632 
0725 
2007 
O656 
2007 
061*6 
2001* 
2006 
2000 
2002 

1655 
l6ll 
1602 
1550 
2002 
15  ^2 
2002 
2001* 
2006 
061*6 
0660 
l6ll 
2000 
061*1* 
11*1*0 


Multiply 
Load 

B  +^/Bz  »1*AC 
Divide  for  >i+1 
Store  in  B 
Load 

Load 

zi  "  zi-l  ■  hi 
Load 

0^i^i)  =  »i+i 

zi 

zi  +  hi+l  ■  zi+l 

Store  in  A 
Load 


zi 
Vzi 


zi+l  -  zi 


^ 


If  e  ^  for  real  part,  0631 
No,  evaluate  f(Zi+i) 
If  t   -q  for  imag  part,  0725 
Build  and 

store  testvord 
Reset  A.  pickup 

in  06*6 
Load 

zi+l 
Load 

An 
Set  exit  in  cpx.  rapy.  routine 

Transfer  into  cpx.  mpy. 

Temp,  store  imag  part 

Set  exit  in  add  routine 

Load  A, 

Transfer  into  add  routine 

Load  imag  part  for  cpx.  mpy. 

Reload  Zi+i 

for  next  term 
Modify  ky   pickup 
If  AQ  not  taken, 

repeat  back 
Testvord  storage 
Constant  for  06**2 
Store  f(si+i) 


TV  -  1  -  7 


LACRANGIAN  INTERPOLATION  METHOD  FOR  HIGH  DEGREE  POLYNOMIALS 


ADD 

INST 

M1 

M* 

M3 

REMARKS 

066l 

51 

1307 

1350 

2004 

Load 

0662 

51 

1527 

1370 

2006 

f(«1)i 

O663 

3^ 

3000 

2100 

1650 

f(zi+1)/f(zi) 

0664 

51 

2000 

2001 

2004 

Load  quotient 

0665 

30 

2002 

2100 

2006 

for  complex 

0666 

56 

2100 

2003 

2007 

conjugate 

0667 

54 

3000 

2100 

1607 

multiplication 

0670 

55 

0676 

2000 

0700 

Ifi 

0671 

55 

2000 

O676 

0673 

|f(zi+1)/f(z<)|2 
is  greater  than 

0672 

5U 

0677 

2001 

0700 

0673 

56 

1550 

1637 

1550 

100;  halve 

067U 

56 

1352 

1637 

1332 

X1+1  and 

0675 

3* 

3000 

2100 

0605 

recompute  2^+1 

0676 

00 

0000 

0000 

0007 

Floating  point 

0677 

00 

5200 

0000 

0650 

100 

0700 

17 

2010 

3000 

0702 

If  SEN  1  up,  print  iteration 

0701 

54 

3000 

2100 

0706 

No,  skip  print 

0702 

21 

1310 

2100 

0001 

Print 

0703 

21 

1551 

2100 

0001 

zi+l 

0704 

21 

1512 

2100 

0001 

real  and 

0705 

21 

1555 

2100 

0001 

imaginary 

0706 

51 

1301 

1342 

1300 

Shift  zj.! 

0707 

51 

1321 

1362 

1320 

to  Zi.g 

0710 

51 

1302 

13^5 

1301 

Shift  zi 

0711 

51 

1322 

1363 

1321 

to  zi-i 

0712 

51 

1310 

1351 

1302 

Shift  zi+1 

0713 

51 

1312 

1353 

1322 

to  Zj 

07l4 

51 

1550 

1371 

1305 

Shift  X1+1 

0715 

51 

1552 

1575 

1325 

tO  \;. 

Shift  ftzi-i) 
to  f  (a^o* 
Shift  ffzi) 

0716 

51 

1506 

15^7 

1305 

0717 

51 

1526 

1567 

1325 

0720 

51 

1507 

1350 

1306 

0721 

51 

1527 

1370 

1326 

to  f (zi_i) 

0722 

51 

1440 

1401 

1307 

Shift  f(zi+i) 

0723 

51 

1442 

1403 

1327 

to  f (Zi) 

0724 

5^ 

5000 

2100 

0361 

Reiterate 

0725 

35 

1512 

1455 

0751 

If  i*ag  portion  significant,  0 

0726 

55 

0751 

1466 

0735 

Initialize  0735 

0727 

52 

0750 

1677 

0737 

Initialize  0757 

0730 

52 

0747 

l64o 

0740 

Initialize  0740 

0731 

31 

0100 

0101 

2000 

Load  Aq 

0732 

25 

2001 

1551 

2005 

}*n  zi 

0735 

55 

2000 

1310 

2002 

0734 

51 

2002 

2003 

2002 

1. 

0755 

51 

0104 

0105 

2000 

An  1  + 

0736 

54 

5000 

2100 

1540 

J     An  zt 

0737 

50 

2000 

2100 

0104 

Put  away 

0740 

50 

2001 

2100 

0105 

terms  of  suppressed  equation 

0741 

55 

0755 

1466 

0755 

Modify  to  suppress 

rv  -  1  -  8 


LAGRANGIAN  INTERPOLATION  METHOD  FOR  HIGH  DEGREE  POLYNOMIALS 


ADD 

n/rr 

M1 

M* 

M3 

0742 

55 

0757 

0264 

0757 

0745 

35 

0740 

0264 

0740 

0744 

3* 

O656 

0755 

0752 

07U5 

36 

1572 

1551 

1572 

0746 

3* 

5000 

2100 

1015 

0747 

00 

0000 

0000 

0105 

0750 

00 

0000 

0000 

0102 

0751 

35 

1310 

1551 

1305 

0752 

50 

1551 

2100 

1344 

0753 

05 

5000 

5000 

1310 

0754 

31 

2002 

2005 

2006 

0755 

30 

2000 

2100 

2004 

0756 

36 

2100 

2001 

2005 

0757 

54 

5000 

2100 

1607 

0760 

31 

2000 

2001 

1504 

0761 

35 

0766 

1475 

2007 

0762 

55 

2007 

0264 

0770 

0763 

32 

1457 

1677 

1002 

0764 

35 

1002 

1551 

2007 

0765 

32 

2007 

1555 

1003 

0766 

31 

0100 

0101 

1300 

0767 

31 

0102 

0103 

1301 

0770 

31 

0114 

0115 

1302 

0771 

35 

1500 

1305 

2000 

0772 

25 

131*1 

1344 

2001 

0773 

31 

2000 

2001 

2004 

0774 

35 

1300 

1504 

2000 

0775 

25 

13^1 

1545 

2001 

0776 

51 

2000 

2001 

2006 

0777 

51 

1301 

1542 

2000 

1000 

51 

1302 

1545 

2002 

1001 

54 

3000 

2100 

1410 

1002 

30 

1300 

2100 

0110 

1003 

30 

154l 

2100 

0111 

1004 

31 

2000 

2001 

1300 

1005 

51 

2002 

2003 

1301 

1006 

35 

0770 

1466 

0770 

1007 

35 

1002 

0264 

1002 

1010 

55 

1005 

0264 

1003 

1011 

54 

0656 

0770 

0770 

1012 

36 

1572 

0264 

1572 

1013 

30 

1551 

1310 

1501 

10l4 

36 

1510 

026l 

2000 

1015 

30 

1551 

2000 

2004 

1016 

55 

2100 

2100 

1502 

1017 

32 

1551 

1166 

1302 

1020 

54 

2004 

2100 

1102 

1021 

5* 

3000 

2100 

1051 

REMARKS 

eqution  to  next 

lover  degree 
If  extraction  not  complete,  0752 
Reduce  N'  by  1 

Transfer  to  conversion  routine 
Constant  for  0750 
Constant  for  0727 
Build  2a  term 

in  quadratic  factor 
Build 

(a2  +  b2) 

term  in 

quadratic 

factor 

Initialize 

0770 
Initialize  1002 
Initialize 

1005 
An  to  TS 

Vl  to  ?S  b 
An-2  t0  TSc 

A^  x  a  to  2004,  2005 


l^   x  (a2  +  b2)  to  2006,  2007 

Load 

Vl 
Add 

Put away  An  of 

suppressed  equation 
Set  up  for 

next  term 
Modify  to  pick  up  Ana=,^ 
Modify  put  away 
Modify  put  away 

Has  extraction  been  completed? 
Yes,  reduce  N'  by  2 
Form  fractional  part   (l) 
Build  shift  control  word 
Form  integer  part 
Clear  sign  storage 
Store  sign 

If  integer  ^  0,  1102 
No,  1051 


fl)  The  remainder  of  Channel  10  and  Channel  11  is  a  binary  floating 
point  to  decimal  conversion  which  would  be  part  of  the  overall  program, 


IV  -  1  -  9 


LAGRANCIAN  INTERPOLATION  METHOD  FOP  HIGH  DEGREE  POLYNOMIALS 


ADD 

INST 

M1 

M* 

M3 

1022 

12 

0000 

0000 

0000 

1023 

10 

0100 

0000 

0000 

102* 

00 

2*00 

0000 

0000 

1025 

00 

1515 

2700 

0000 

1026 

10 

0000 

0100 

0000 

1027 

00 

3115 

1515 

1500 

1030 

30 

2005 

2100 

1300 

1031 

30 

1301 

2100 

2000 

1032 

30 

1167 

2100 

2007 

1033 

26 

2000 

0257 

2000 

103* 

27 

2007 

0256 

2007 

1035 

32 

2001 

025* 

2007 

IO36 

37 

2007 

2100 

1033 

1037 

50 

2007 

2100 

1301 

10*0 

32 

1016 

1022 

1301 

10*1 

3* 

II63 

1067 

10*3 

10*2 

5* 

3000 

2100 

1170 

10*3 

3* 

1331 

200* 

1057 

10** 

27 

1101 

2100 

2007 

10*5 

32 

1300 

0260 

2100 

10*6 

3* 

2000 

2100 

1055 

10*7 

21 

102* 

1023 

0001 

1050 

30 

1300 

0256 

1300 

1051 

30 

2007 

0256 

2007 

1052 

3* 

3000 

2100 

10*5 

1053 

32 

1302 

1022 

1300 

105* 

21 

1500 

2007 

0001 

1055 

21 

1025 

1026 

0001 

1056 

3* 

3000 

2100 

1066 

1057 

21 

1027 

11*0 

0001 

1060 

32 

1302 

1166 

2100 

106l 

27 

2000 

0262 

2000 

1062 

3* 

2000 

2100 

IO65 

IO63 

21 

1157 

1026 

0001 

106* 

3* 

3000 

2100 

1066 

IO65 

21 

ll60 

1026 

0001 

1066 

21 

1301 

1101 

0001 

IO67 

55 

1455 

1312 

107* 

1070 

21 

1161 

II65 

0002 

1071 

35 

1163 

2100 

1067 

1072 

31 

1312 

1553 

1310 

1073 

3* 

3000 

2100 

1013 

107U 

21 

115* 

1023 

0001 

1075 

3* 

1372 

2100 

0310 

1076 

22 

0000 

0000 

0000 

1077 

35 

116* 

2100 

1067 

1100 

3* 

3000 

2100 

107* 

1101 

01 

0000 

0000 

0001 

REMARKS 


Suppress  sign  bits 

Print  control 

Alphabetic  print 

Alphabetic  print 

Print  control 

Alphabetic  print 

Reset  command 

Load  fraction 

Load  testword 

Multiply  for  most  sig.  digit 

Shift  testword 

Extract  digit  into  TV  storage 

Repeat  loop  if  o'flo 

Putaway  converted  fraction 

Suppress  sign 

Print  editing 

commands 
Print  editing 

commands 
Extract  most  sig.  dec.  digit 
If  2000  *  0,  1053 
If  not  print  space 
Shift  word  *  left 

and  change  print  control 
Return  to  test  next  digit 
Enter  sign 
Print  integer 
Print  decimal  point 
Skip  to  1066 
Print  editing  command 
Build  sign  print 

when  only  fractional 

part  involved 
Print  4. 
Skip  to  1066 
Print  -. 

Print  fractional  part 
If  imag  part  insig,  107* 
No,  print  +  j 
Set  skip  command  in  IO67 
Set  up  convert 

imaginary  part 
Print  carriage  return 
If  N'  ^  0,  next  root 
If  IT  =0,  halt 
Restore  1067 

Return  to  "test  for  end* 
Constant  for  10** 


IV  -  1  -  10 


LAGRANGIAN  INTERPOLATION  METHOD  FOR  HIGH  DEGREE  POLYNOMIALS 


ADD 

INST 

M1 

M2 

M3 

REMARKS 

1102 

35 

2100 

2100 

2005 

High 

1103 

31 

1113 

2004 

2000 

integer 

1104 

32 

2000 

1114 

2100 

conversion 

1105 

30 

2000 

1115 

2000 

routine 

1106 

35 

1116 

2000 

1107 

1107 

35 

1126 

2100 

1110 

1110 

23 

2004 

1140 

2000 

1111 

26 

2000 

1147 

2000 

1112 

3l* 

3000 

2100 

1142 

1113 

02 

0000 

0000 

0002 

1114 

00 

0000 

0000 

0174 

1115 

00 

0000 

0000 

0026 

1116 

35 

1115 

2100 

1110 

1117 

23 

2004 

1127 

2000 

1120 

23 

2004 

1130 

2000 

1121 

23 

2004 

1131 

2000 

1122 

23 

2004 

1132 

2000 

1123 

23 

2004 

1133 

2000 

1124 

23 

2004 

1134 

2000 

1125 

23 

2004 

1136 

2000 

1126 

23 

2004 

1137 

2000 

1127 

00 

0167 

1531 

1777 

1130 

00 

0013 

3274 

0777 

1131 

00 

0001 

1422 

2377 

1132 

00 

0000 

0750 

2177 

1133 

00 

0000 

0060 

6477 

1134 

00 

0000 

0004 

7037 

1135 

34 

3000 

2100 

1151 

1136 

00 

0000 

0000 

0307 

1137 

00 

0000 

0000 

0023 

ll4o 

10 

0000 

0000 

0100 

n4i 

26 

2000 

1155 

2000 

11U2 

30 

2005 

1150 

2005 

11U3 

35 

1110 

1156 

1110 

n44 

35 

2001 

2005 

2005 

11 45 

33 

1110 

1153 

ll4l 

111*6 

34 

3000 

2100 

1030 

11U7 

00 

0000 

0000 

0024 

1150 

00 

0000 

0000 

0004 

1151 

35 

2100 

2100 

2005 

1152 

34 

3000 

2100 

1146 

1153 

23 

2004 

1140 

2000 

1154 

00 

3400 

0000 

0000 

1155 

00 

0000 

OCOO 

0012 

1156 

00 

0000 

0001 

0000 

IV  -  1  -  11 


LAGRANCIAN  INTERPOLATION  METHOD  FOR  HIGH  DEGREE  POLYNOMIALS 


ADD 

INST 

M1 

M* 

M3 

1157 

00 

2352 

276U 

61*00 

ll60 

00 

2252 

2761* 

61*00 

ll6l 

00 

2315 

3621 

3261* 

1162 

00 

k26k 

61*64 

3232 

II63 

3U 

3000 

2100 

1077 

116U 

33 

1U55 

1312 

107** 

II65 

10 

0000 

0000 

0000 

1166 

02 

0000 

0000 

0000 

II67 

00 

0U21 

01*21 

0U20 

1170 

32 

1016 

1022 

1302 

1171 

34 

200U 

2100 

101*1* 

1172 

21 

Ulh 

1175 

0001 

1175 

5k 

3000 

2100 

1066 

117U 

00 

5227 

0000 

0000 

1175 

10 

0001 

0000 

0000 

1176  ■ 

-  1177 

REMARKS 


Print  configuration  for  IO63 
Print  configuration  for  IO65 
Print  configuration  for  1070 
Print  configuration  for  1070 
Dummy  reset  for  1071 
Dummy  reset  for  107** 
Print  control  -  alphabetic 
Sign  extractor 
Testword  for  1072 
Suppress  integer  sign  when 

printing  imaginary  portion 

or  print  0. 

and  return  to  1066 
Print  configuration  for  1172 
Print  control  word 
Not  required 


Channel  12  is  used  for  a  memory  check- sum  routine  in  the  root -solving 
routine  and  is  available  for  other  use  when  this  program  is  incorporated 
in  the  overall  synthesis  program. 

Channel  13  is  used  for  temporary  storage  and  for  a  few  constants  which 
are  listed  below.  NOTE:  This  program  is  coded  for  minimum  access  only, 
i.e.,  adjacent  cells  in  the  main  memory  are  numbered  00,  1*1,  02,  1*3, 
0k,   1*5,  etc. 


1331,  1357  -  floating  point  "1",  must  be  addressed  Individually 
1333,  1335  -  floating  point  "0",  must  be  addressed  individually 


1372  -  N*  -  the  degree  of  the  reduced  equation 

1371*   00    3000    3000    0200       Testword  for  0305 

1376   00    0000    0000    0010       Increment  for  O303,  O30U 


Channels  lU,  15,  and  16  contain  complex  floating  point  arithmetic 
routines,  i.e., 


Complex  Addition 
Complex  Subtraction 
Simple  Multiplication 
Simple  Division 
Simple  Subtraction 
Simple  Addition 
Square  Root 

Complex  Multiplication 
Complex  Division 


11*10 
H*l6 
1500 
1507 
1516 
15U0 
1551 
1607 
1650 


IV  -  1  -  12 


LAGRANGIAN  INTERPOLATION  METHOD  FOR  HIGH  DEGREE  POLYNOMIALS 


ADD 

INST 

M1 

M* 

M3 

REMARKS 

1400 

-  1407 

Used 

for  temporary  storage 

of  operands 

i4io 

30 

3000 

1505 

i4oo 

Prepare 

l4ll 

32 

1400 

l64o 

1433 

exit 

l4l2 

04 

3000 

3000 

1400 

Store  operands 

1U13 

32 

1U72 

1535 

1424 

Set  1424 

l4l4 

32 

1U71 

1677 

1430 

and  1430  for  addition 

l4l5 

34 

3000 

2100 

1423 

Transfer  to  1423 

l4l6 

30 

3000 

1603 

1536 

Prepare 

1417 

32 

1536 

1515 

1433 

exit 

1420 

04 

3000 

3000 

1400 

Store  operands 

1421 

32 

IU65 

I636 

1424 

Set  1424 

1422 

32 

1U67 

1636 

1430 

and  1430  for  subtraction 

1423 

31 

2004 

2005 

2002 

Reposition  real  operands 

1424 

?k 

3000 

2100 

1516 

Transfer  for  ♦  reals 

1U25 

51 

2000 

2001 

1400 

Store  sum  of  reals 

1426 

31 

l402 

1443 

2000 

Load  imaginary 

1U27 

31 

l4o6 

1447 

2002 

operands 

1U30 

34 

3000 

2100 

1516 

Transfer  for  +  imaginaries 

IU31 

31 

2000 

2001 

2002 

Reposition  Imaginary  sum 

1^52 

31 

1400 

i44i 

2000 

Load  real  sum 

1433 

3k 

3000 

2100 

0624 

Exit 

1434  ■ 

-  1U36 

Temporary  storage 

1437 

00 

0000 

0000 

0100 

Decrementer  for  0306,  0307 

1440  ■ 

■  1453 

Temporary  storage 

1U5U 

02 

0000 

0000 

0044 

1455  « 

.  1W63 

Temporary  storage 

1464 

00 

0000 

7777 

7777 

M2,  M3  extractor  for  0333 

1465 

00 

0000 

0000 

1516 

Constant  for  1421 

1466 

00 

0002 

0002 

0000 

Constant  for  0323 

1467 

00 

0000 

0000 

1516 

Constant  for  1422 

1U70 

31 

0102 

0103 

2000 

Dummy  reset  for  0326 

1U71 

00 

0000 

0000 

1540 

Constant  for  l4l4 

1U72 

00 

0000 

0000 

1540 

Constant  for  l4l3 

1473 

00 

ooo4 

ooo4 

0000 

Constant  for  03 32,  O336 

1U7U 

00 

0002 

0000 

0000 

Incrementer  (various ) 

1U75  - 

■  1477 

Temporary  storage 

1500 

30 

3000 

1505 

2006 

Prepare 

1301 

32 

2006 

1677 

1550 

exit 

1502 

35 

2000 

2004 

2000 

Add  exponents 

1503 

25 

2001 

2005 

2001 

Multiply  magnitudes 

150U 

5k 

3000 

2100 

15^7 

Transfer  to  exit 

1505 

02 

0000 

0000 

0030 

Shift  control  for  1500 

1506 

00 

0000 

0000 

3777 

M3  extractor 

1507 

30 

3000 

1676 

2006 

Prepare 

1510 

32 

2006 

1677 

1550 

exit 

1511 

36 

2000 

2004 

2000 

Subtract  exponents 

1512 

23 

2001 

2005 

2001 

Divide  magnitudes 

1515 

37 

2001 

3000 

1530 

If  overflow,  1530 

15l4 

5k 

3000 

2100 

15^7 

Transfer  to  exot 

1515 

00 

0000 

0000 

7777 

M3  extractor 

IV  -  1  -  15 


LAGRANOIAN  INTERPOLATION  METHOD  FOR  HIGH  DEGREE  POLYNOMIALS 


ADD 

?r,?T 

M1 

Mz 

M3 

1516 

30 

3000 

1505 

2006 

1517 

32 

2006 

1677 

1550 

1520 

36 

2100 

2003 

2003 

1521 

3^ 

3000 

2100 

15^2 

1522 

36 

2000 

2002 

2006 

1523 

30 

2001 

2006 

2007 

152U 

50 

2002 

2100 

2000 

1525 

35 

2003 

2007 

2001 

1526 

57 

2001 

3000 

1550 

1527 

54 

3000 

2100 

15^7 

1530 

27 

2001 

l6kk 

2001 

1531 

30 

2001 

l6kk 

2001 

1532 

27 

2001 

I6k6 

2001 

1533 

35 

2000 

l6k6 

2000 

153^ 

3^ 

3000 

2100 

I5U7 

1535 

00 

0000 

0000 

7777 

1536 

00 

0000 

0000 

062U 

1537 

02 

0000 

0000 

0001 

15*K) 

30 

3000 

1505 

2006 

15^1 

32 

2006 

1515 

1550 

15U2 

33 

2002 

2000 

1522 

15^3 

36 

2002 

2000 

2006 

15  M* 

30 

2003 

2006 

2007 

15^5 

35 

2001 

2007 

2001 

15^6 

37 

2001 

3000 

1530 

15^7 

31 

2000 

2001 

2000 

1550 

3»* 

3000 

2100 

•  •  •  • 

1551 

30 

3000 

1676 

2006 

1552 

32 

2006 

1535 

1550 

1553 

35 

2000 

026l 

2000 

1554 

32 

2000 

1331 

210U 

1555 

30 

2000 

l6kk 

2006 

1556 

3k 

200U 

2100 

1561 

1557 

30 

2001 

l6kk 

2007 

1560 

3k 

3000 

2100 

1562 

1561 

30 

2001 

2100 

2007 

1562 

50 

1575 

2100 

2002 

1563 

23 

2007 

2002 

2001 

156U 

36 

2001 

2002 

2001 

1565 

30 

2001 

l6kh 

2001 

1566 

3^ 

2001 

2100 

1570 

1567 

3k 

3000 

2100 

1572 

1570 

35 

2002 

2001 

2002 

1571 

3k 

3000 

2100 

1563 

1572 

25 

2002 

1577 

2001 

1573 

36 

2006 

1576 

2000 

1574 

3k 

3000 

2100 

15^7 

REMARKS 


>nd 


Prepare 

exit 
Change  sign  2"u  operand 
Transfer  into  add 
Form  shift  control 
Shift  Is*  operand 
Shift  2   exponent  to  Ans.  cell 
Add  operands 
If  overflow,  1530 
Transfer  to  exit 
Normal 
correct 

overflow 
procedure 
Transfer  to  exit 
M3  extractor 
Prepare  exit  storage 
Shift  control  for  1555 
Prepare 

exit 
If  2002  ^  2000,  1522 
No,  form  shift  control 
Shift  2nd  operand 
Add  operands 
If  overflow,  1530 
Normalize 
Exit 
Prepare  exit 

from  square  root  routine 
Multiply  by  244 

Extract  last  bit  of  exp  into  2004 
Take  square  root  of  exp 
If  exp  odd,  magnitude  has  been 

divided  by  2,  by  1555 
If  exp  even,  divide  magnitude 

by  2 
Skip  to  iteration 
Load  magnitude  for  iteration 
Load  a.* 
(x/2)/&1 
fx/2a,  )  -  at 

l/2(x/2ai)  -  b.±   =  ai+i  -  at 
Has  process  converged? 
Yes,  1572;  no,  1570 

ai+l  r^Pl810®8  ai 
and  reiterate 

Multilpy  by  42/2 

Divide  by  221 

Transfer  to  exit 


IV  -  1  -  lU 


LAGRANCIAN  INTERPOLATION  METHOD  FOR  HIGH  DEGREE  POLYNOMIALS 


ADD 

INST 

Mx 

M2 

M3 

REMARKS 

1575 

00 

7777 

7777 

7777 

Constant  for  1562 

1576 

00 

0000 

0000 

0021 

Constant  for  1573 

1577 

00 

5520 

2363 

1500 

Constant  for  1572 

1600 

-  1602 

Used 

for  temporary  storage 

of  operands 

1603 

02 

0000 

0000 

0030 

Shift  control  for  1607 

l6o4 

-  1606 

Used 

for  temporary  storage 

of  operands 

1607 

30 

3000 

1603 

1601 

Prepare 

1610 

32 

1601 

1636 

1635 

exit 

1611 

04 

3000 

3000 

l600 

Store  operands 

1612 

34 

5000 

2100 

1500 

Multiply  reals 

1613 

30 

2000 

2100 

1601 

and  store  results 

l6l4 

30 

2001 

2100 

1605 

in  l601,  1605 

1615 

55 

1602 

1606 

2000 

Multiply  imaginaries 

1616 

25 

1643 

1647 

2001 

to  obtain  negative  real 

1617 

31 

2000 

2001 

2002 

and  shift  for  subtraction 

1620 

31 

1601 

l605 

2000 

from  previous  real  product 

1621 

3^ 

3000 

2100 

1516 

Transfer  to  subtract 

1622 

30 

2000 

2100 

1601 

Putaway  real  product 

1623 

30 

2001 

2100 

1605 

in  1601,  1605 

1624 

35 

1602 

1604 

2000 

Multiply  for 

1625 

25 

1643 

1645 

2001 

first  cross-product  term 

1626 

31 

2000 

2001 

2002 

and  hold  in  2002,  2003 

1627 

35 

1600 

1606 

2000 

Multiply  second 

1630 

25 

l64l 

1647 

2001 

cross-product  term 

1631 

31 

2000 

2001 

2000 

and  prepare  for  addition 

1632 

3U 

3000 

2100 

1540 

Add  for  imaginary  cross-products 

1633 

31 

2000 

2001 

2002 

Put  results  in2002,  2003 

1634 

31 

1601 

1605 

2000 

Insert  real  product 

1635 

54 

3000 

2100 

0627 

Exit 

1636 

00 

0000 

0000 

7777 

M3  extractor  for  1610 

1637 

00 

0000 

0000 

0001 

Shifter  for  1656 

l64o 

00 

0000 

0000 

7777 

M3  extractor 

l64i 

Used 

for  temporary  1 

storage 

1642 

00 

0000 

0000 

0001 

Shift  control  word 

1643 

Used 

for  temporary  : 

storage 

1644 

02 

0000 

0000 

0001 

Shift  control  word 

1645 

Used 

for  temporary  storage 

1646 

00 

0000 

0000 

0001 

Shift  Control  word 

16U7 

Used 

for  temporary  storage 

1650 

30 

3000 

1676 

1601 

Prepare      (A  +JB) 

1651 

32 

1601 

1677 

1635 

exit      (C  +JD) 

1652 

04 

3000 

3000 

1600 

Store  operands 

1653 

50 

2004 

1637 

2000 

''orm 

1654 

25 

2005 

1645 

2001 

C2 

1655 

31 

2000 

2001 

2000 

and  normalize 

TV  -  1  -  15 


LAGRANOIAN  INTERPOLATION  METHOD  FOR  HIGH  DEGREE  POLYNOMIALS 


ADD   INST 


MJ 


M2 


M^ 


REMARKS 


16^6   30 
1657   25 


1660 
1661 
1662 
I663 
l664 
1665 


1670 
1671 
1672 
1673 
I67U 

1675 
1676 

1677 


31 
34 
51 

31 
34 

30 


1666  30 

1667  31 


34 
30 
36 
31 
31 
34 
02 
00 


2006 
2007 
2002 
3000 
2000 
l604 
3000 
2000 
2001 
1606 
3000 
2000 
2100 
l604 
1600 
3000 
0000 
0000 


1637 
I6U7 
2003 
2100 
2001 
16U5 
2100 
2100 
2100 
16V7 
2100 
2100 
2001 

1645 
l64l 

2100 
0000 
0000 


2002 
2003 
2002 
15^0 
2004 
2000 

1507 
1604 

1645 
2000 

1507 
1606 

1647 
2004 
2000 
1612 

0030 

7777 


Form 

D2 

end  normalize 
(C2  +  D2) 

Shift  to  2004,  2005 
Load  C 

C/  (C2  +  D2) 
Store  in 

l604,  1645 
Load  D 
D/  (C2  +  D2) 
Store  in  1606, 

1647  vith  negative  sign 
Load  conjugate 

quotient  and 

transfer  to  multiply 
Shift  control  word 
M3  extractor 


FINIS 


IV  -  1  -  16 


FLOW  CHART  OF  LAG RANG UN   INTERPOLATION  METHOD 
FOR  SOLUTION  OF  ALGEBRAIC  EQUATIONS   OF  HIGH  DEGREE 


(Decimal  Coefficient   In  CH  00   )  (  Degree  of  Pn(z)      ) 


Conv.    to  Bin. Fit.   Pt.    «£ 


Coeffs    to  CH  01 

N— =>  N' 


-1  +  JO 
1  +  JO 
0  +  JO 

4  +  j° 


f<Zl-2> 

f<Vl> 
f(zt) 


Xt 


1+1 


1-1,2,3. .n/2 


21 
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APPENDIX  V 
USE  OF  THE  EXISTING  PROGRAMS  ON  THE  CRC-102A 

All  existing  programs  have  been  coded  for  minimum-access  word  channel 
and  will  not  work  in  the  sequential  channels. 

1.  Dirichlet  Power  Series  Expansion 

Since  the  synthetic  division  and  conversion  program  has  not  been 
written  as  yet,  if  f  (t)  and  f.(t)  are  given,  \h\        must  be  computed  by 
hand.   If  h(t)  is  given,  it  must  be  converted  to  ijh.3   by  multiplying  by 
At.   As  the  program  is  now  located,  storage  is  available  for  100  (octal) 
ordinates  in  floating  point  form  if  the  decimal  to  floating  point  con- 
version in  Channel  2  is  to  be  used.   If  the  conversion  is  accomplished 
prior  to  entry  of  data,  150  (octal)  ordinates  may  be  used. 

T 
S  ■  2  -—  a  2NT      where  S  ■  total  cells  needed 
At 

T  s  number  of  time  units  problem  is 
scaled  over  (note  command  0423) 

At  a»  sampling  interval 

N  ■  samples  per  time  unit 

The  program  may  be  relocated  to  accommodate  more  sampled  ordinates. 

The  binary  floating  point  numbers  representing  the  ordinates  of  h(t) 
are  entered  in  cells  0000  -  0117,  exponents  in  even  cells,  magnitudes  in 
odd;  this  allows  40  ordinates  to  be  used. 

The  test  switches  determine  the  amount  of  interpolation  performed. 
The  output  will  be  H*(s)  in  decimal  form  for  as  many  terms  requested. 
The  number  of  terms  is  set  by  entering  this  datum  as  an  octal  number  in 
J,  cell  0540. 

2.  Matrix  Inversion 

The  matrix  inversion  program  takes  the  ascending  power  H*(s)  terms  in 
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packed  floating  point  form  (a  conversion  routine  is  available)  in  the  word 
structure  shown  on  page  III -5  and  computes  the  lossless  system  function, 
Z'(p).   The  terms  are  entered  in  Channel  0,0000,  0001,  0002,  etc.   Cells 
0076  and  0077  are  reserved  for  N,  the  order  of  the  matrix  which  equals  n, 
the  number  of  poles  desired  or  the  degree  of  the  denominator  and  n  ■♦•  r, 
the  degree  of  the  numerator  respectively.   For  the  example  on  page  21, 
N  »  n  «  h  (0076) 

n  +  r=n-2-2   (0077) 

The  commands  0300  through  0320  contain  instructions  which  link  the  program 
to  a  packed  binary  floating  point  to  decimal  conversion  routine  located  in 
1571  -  1677. 

If  these  commands  are  included,  the  program  will  yield  the  decimal 
coefficients  of  Z'p  in  the  following  form: 

pQ  +  p  8  +  P2S2 


q0  +  qlS  +  q2S  *  q3S   +  q4S 


qn  will  be  equal  to  1,  the  results  must  be  normalized  to  the  following 
form  for  entry  in  the  root-solving  routine. 

2 
(M)    s   +  p|s  +  P(J 

~4     "3      2 

s   +  q^s   +  q^s   +  qjs  +  q^ 

3.   Lagrangian  Interpolation  Method 

The  procedure  for  the  use  of  this  program  is  clearly  described  in 
the  description  of  sub-routines  for  the  Computation  Center  (Routine  2.1). 

It  is  recommended  that  the  Channel  2  conversion  routine  of  this  pro- 
gram be  replaced  to  provide  the  normalizing  procedure  required  by  2.  above 
Channels  11  and  12  are  available  for  the  storage  of  the  factored  forms  of 
the  numerator  and  denominator. 
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